diff options
Diffstat (limited to 'Examples/Modules')
23 files changed, 632 insertions, 506 deletions
diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py index 3dcf40ac8..f28272897 100755 --- a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py +++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py @@ -19,8 +19,8 @@ import read_raw_data yt.funcs.mylog.setLevel(0) # Read data from back-transformed diagnostics -snapshot = './lab_frame_data/snapshot00001' -header = './lab_frame_data/Header' +snapshot = './lab_frame_data/snapshots/snapshot00001' +header = './lab_frame_data/snapshots/Header' allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) z = np.mean( read_raw_data.get_particle_field(snapshot, 'beam', 'z') ) w = np.std ( read_raw_data.get_particle_field(snapshot, 'beam', 'x') ) diff --git a/Examples/Modules/RigidInjection/inputs.BoostedFrame b/Examples/Modules/RigidInjection/inputs.BoostedFrame index c7a60f14f..456df363d 100644 --- a/Examples/Modules/RigidInjection/inputs.BoostedFrame +++ b/Examples/Modules/RigidInjection/inputs.BoostedFrame @@ -3,7 +3,7 @@ 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.do_back_transformed_diagnostics = 1 warpx.num_snapshots_lab = 2 warpx.dt_snapshots_lab = 1.8679589331096515e-13 diff --git a/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py b/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py new file mode 100755 index 000000000..3b8e7aa76 --- /dev/null +++ b/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py @@ -0,0 +1,36 @@ +#! /usr/bin/env python + +''' +Analysis script of a WarpX simulation in a boosted frame. + +The simulation runs in a boosted frame, and the analysis is done in the lab +frame, i.e., on the back-transformed diagnostics for the full 3D simulation and +an x-z slice at y=y_center. The field-data, Ez, along z, at (x_center,y_center,:) is compared +between the full back-transformed diagnostic and the reduced diagnostic (i.e., x-z slice) . +''' + +import numpy as np +import read_raw_data + +# Read data from back-transformed diagnostics of entire domain +snapshot = './lab_frame_data/snapshots/snapshot00002' +header = './lab_frame_data/snapshots/Header' +allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) +F = allrd['Ez'] +print("F.shape ", F.shape) +F_1D = np.squeeze(F[F.shape[0]//2,F.shape[1]//2,:]) + + +# Read data from reduced back-transformed diagnostics (i.e. slice) +snapshot_slice = './lab_frame_data/slices/slice00002' +header_slice = './lab_frame_data/slices/Header' +allrd, info = read_raw_data.read_lab_snapshot(snapshot_slice, header_slice) +Fs = allrd['Ez'] +print("Fs.shape", Fs.shape) +Fs_1D = np.squeeze(Fs[Fs.shape[0]//2,1,:]) + +error = np.max(np.abs(Fs_1D - F_1D)) / np.max(np.abs(F_1D)) + +# Print error and assert small error +print("relative error: " + str(error)) +assert( error < 1E-15 ) diff --git a/Examples/Modules/boosted_diags/inputs.2d b/Examples/Modules/boosted_diags/inputs.2d deleted file mode 100644 index 6afe6977d..000000000 --- a/Examples/Modules/boosted_diags/inputs.2d +++ /dev/null @@ -1,95 +0,0 @@ -# Maximum number of time steps -max_step = 260 - -# number of grid points -amr.n_cell = 64 64 512 - -# Maximum allowable size of each subdomain in the problem domain; -# this is used to decompose the domain for parallel calculations. - -amr.max_grid_size = 32 - -# Maximum level in hierarchy (for now must be 0, i.e., one level in total) -amr.max_level = 0 -amr.plot_int = 10 # How often to write plotfiles. "<= 0" means no plotfiles. -amr.check_int = 10 - -# Geometry -geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? -geometry.prob_lo = -150.e-6 -150.e-6 -0.6e-3 # physical domain -geometry.prob_hi = 150.e-6 150.e-6 0. - -# Verbosity -warpx.verbose = 1 - -# Algorithms -algo.current_deposition = direct - -# Numerics -interpolation.nox = 3 -interpolation.noy = 3 -interpolation.noz = 3 -warpx.use_filter = 1 -warpx.cfl = 1.0 -warpx.do_pml = 0 - -# Moving window -warpx.do_moving_window = 1 -warpx.moving_window_dir = z -warpx.moving_window_v = 1.0 # in units of the speed of light - -# Boosted frame -warpx.gamma_boost = 15. -warpx.boost_direction = z - -# Diagnostics -warpx.do_boosted_frame_diagnostic = 1 -warpx.num_snapshots_lab = 20 -warpx.dt_snapshots_lab = 7.0e-14 - -# Species -particles.nspecies = 2 -particles.species_names = electrons ions - -electrons.charge = -q_e -electrons.mass = m_e -electrons.injection_style = "NUniformPerCell" -electrons.xmin = -150.e-6 -electrons.xmax = 150.e-6 -electrons.ymin = -150.e-6 -electrons.ymax = 150.e-6 -electrons.zmin = 0.e-6 -electrons.num_particles_per_cell_each_dim = 1 1 2 -electrons.profile = constant -electrons.density = 1. -electrons.momentum_distribution_type = "constant" -electrons.do_continuous_injection = 1 - -ions.charge = q_e -ions.mass = m_p -ions.injection_style = "NUniformPerCell" -ions.xmin = -150.e-6 -ions.xmax = 150.e-6 -ions.ymin = -150.e-6 -ions.ymax = 150.e-6 -ions.zmin = 0.e-6 -ions.num_particles_per_cell_each_dim = 1 1 2 -ions.profile = constant -ions.density = 1. -ions.momentum_distribution_type = "constant" -ions.do_continuous_injection = 1 - -# Laser -lasers.nlasers = 1 -lasers.names = laser1 -laser1.profile = Gaussian -laser1.position = 0. 0. -1.e-6 # This point is on the laser plane -laser1.direction = 0. 0. 1. # The plane normal direction -laser1.polarization = 1. 0. 0. # The main polarization vector -laser1.e_max = 8.e12 # Maximum amplitude of the laser field (in V/m) -laser1.profile_waist = 5.e-5 # The waist of the laser (in meters) -laser1.profile_duration = 16.7e-15 # The duration of the laser (in seconds) -laser1.profile_t_peak = 33.4e-15 # The time at which the laser reaches its peak (in seconds) -laser1.profile_focal_distance = 0.e-6 # Focal distance from the antenna (in meters) -laser1.wavelength = 0.8e-6 # The wavelength of the laser (in meters) diff --git a/Examples/Modules/boosted_diags/inputs.3d b/Examples/Modules/boosted_diags/inputs.3d deleted file mode 100644 index 528eb6cd9..000000000 --- a/Examples/Modules/boosted_diags/inputs.3d +++ /dev/null @@ -1,95 +0,0 @@ -# Maximum number of time steps -max_step = 260 - -# number of grid points -amr.n_cell = 64 64 512 - -# Maximum allowable size of each subdomain in the problem domain; -# this is used to decompose the domain for parallel calculations. - -amr.max_grid_size = 32 - -# Maximum level in hierarchy (for now must be 0, i.e., one level in total) -amr.max_level = 0 -amr.plot_int = 10 # How often to write plotfiles. "<= 0" means no plotfiles. -amr.check_int = 10 - -# Geometry -geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 1 1 0 # Is periodic? -geometry.prob_lo = -150.e-6 -150.e-6 -0.6e-3 # physical domain -geometry.prob_hi = 150.e-6 150.e-6 0. - -# Verbosity -warpx.verbose = 1 - -# Algorithms -algo.current_deposition = direct - -# Numerics -interpolation.nox = 3 -interpolation.noy = 3 -interpolation.noz = 3 -warpx.use_filter = 1 -warpx.cfl = 1.0 -warpx.do_pml = 0 - -# Moving window -warpx.do_moving_window = 1 -warpx.moving_window_dir = z -warpx.moving_window_v = 1.0 # in units of the speed of light - -# Boosted frame -warpx.gamma_boost = 15. -warpx.boost_direction = z - -# Diagnostics -warpx.do_boosted_frame_diagnostic = 1 -warpx.num_snapshots_lab = 20; -warpx.dt_snapshots_lab = 7.0e-14; - -# Species -particles.nspecies = 2 -particles.species_names = electrons ions - -electrons.charge = -q_e -electrons.mass = m_e -electrons.injection_style = "NUniformPerCell" -electrons.xmin = -150.e-6 -electrons.xmax = 150.e-6 -electrons.ymin = -150.e-6 -electrons.ymax = 150.e-6 -electrons.zmin = 0.e-6 -electrons.num_particles_per_cell_each_dim = 1 1 2 -electrons.profile = constant -electrons.density = 1. -electrons.momentum_distribution_type = "constant" -electrons.do_continuous_injection = 1 - -ions.charge = q_e -ions.mass = m_p -ions.injection_style = "NUniformPerCell" -ions.xmin = -150.e-6 -ions.xmax = 150.e-6 -ions.ymin = -150.e-6 -ions.ymax = 150.e-6 -ions.zmin = 0.e-6 -ions.num_particles_per_cell_each_dim = 1 1 2 -ions.profile = constant -ions.density = 1. -ions.momentum_distribution_type = "constant" -ions.do_continuous_injection = 1 - -# Laser -lasers.nlasers = 1 -lasers.names = laser1 -laser1.profile = Gaussian -laser1.position = 0. 0. -1.e-6 # This point is on the laser plane -laser1.direction = 0. 0. 1. # The plane normal direction -laser1.polarization = 1. 0. 0. # The main polarization vector -laser1.e_max = 8.e12 # Maximum amplitude of the laser field (in V/m) -laser1.profile_waist = 5.e-5 # The waist of the laser (in meters) -laser1.profile_duration = 16.7e-15 # The duration of the laser (in seconds) -laser1.profile_t_peak = 33.4e-15 # The time at which the laser reaches its peak (in seconds) -laser1.profile_focal_distance = 0.e-6 # Focal distance from the antenna (in meters) -laser1.wavelength = 0.8e-6 # The wavelength of the laser (in meters) diff --git a/Examples/Modules/boosted_diags/inputs.3d.slice b/Examples/Modules/boosted_diags/inputs.3d.slice new file mode 100644 index 000000000..408b18931 --- /dev/null +++ b/Examples/Modules/boosted_diags/inputs.3d.slice @@ -0,0 +1,112 @@ +warpx.zmax_plasma_to_compute_max_step = 0.0031 + +amr.n_cell = 32 32 64 +amr.max_grid_size = 64 +amr.blocking_factor = 32 +amr.max_level = 0 +amr.plot_int = 10000 + +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 1 0 # Is periodic? +geometry.prob_lo = -128.e-6 -128.e-6 -40.e-6 +geometry.prob_hi = 128.e-6 128.e-6 0.96e-6 + +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = energy-conserving +algo.particle_pusher = vay +algo.maxwell_fdtd_solver = ckc +interpolation.nox = 3 +interpolation.noy = 3 +interpolation.noz = 3 +warpx.use_filter = 1 +warpx.cfl = 1. +warpx.do_pml = 0 + +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 + +warpx.gamma_boost = 10. +warpx.boost_direction = z +warpx.do_back_transformed_diagnostics = 1 +warpx.num_snapshots_lab = 4 +warpx.dz_snapshots_lab = 0.001 +warpx.boosted_frame_diag_fields= Ex Ey Ez By rho + +particles.nspecies = 3 +particles.species_names = electrons ions beam +particles.use_fdtd_nci_corr = 1 + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = NUniformPerCell +electrons.num_particles_per_cell_each_dim = 1 1 1 +electrons.momentum_distribution_type = "gaussian" +electrons.xmin = -120.e-6 +electrons.xmax = 120.e-6 +electrons.ymin = -120.e-6 +electrons.ymax = 120.e-6 +electrons.zmin = 0. +electrons.zmax = .003 +electrons.profile = constant +electrons.density = 3.5e24 +electrons.do_continuous_injection = 1 +electrons.do_back_transformed_diagnostics = 1 + +ions.charge = q_e +ions.mass = m_p +ions.injection_style = NUniformPerCell +ions.num_particles_per_cell_each_dim = 1 1 1 +ions.momentum_distribution_type = "gaussian" +ions.xmin = -120.e-6 +ions.xmax = 120.e-6 +ions.ymin = -120.e-6 +ions.ymax = 120.e-6 +ions.zmin = 0. +ions.zmax = .003 +ions.profile = constant +ions.density = 3.5e24 +ions.do_continuous_injection = 1 +ions.do_back_transformed_diagnostics = 1 + +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 = .2e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -20.e-6 +beam.npart = 1000 +beam.q_tot = -1.e-14 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 200000. +beam.ux_th = .2 +beam.uy_th = .2 +beam.uz_th = 20. + +lasers.nlasers = 1 +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0. 0. -0.1e-6 # This point is on the laser plane +laser1.direction = 0. 0. 1. # The plane normal direction +laser1.polarization = 0. 1. 0. # The main polarization vector +laser1.e_max = 2.e12 # Maximum amplitude of the laser field (in V/m) +laser1.profile_waist = 45.e-6 # The waist of the laser (in meters) +laser1.profile_duration = 20.e-15 # The duration of the laser (in seconds) +laser1.profile_t_peak = 40.e-15 # The time at which the laser reaches its peak (in seconds) +laser1.profile_focal_distance = 0.5e-3 # Focal distance from the antenna (in meters) +laser1.wavelength = 0.81e-6 # The wavelength of the laser (in meters) + +slice.dom_lo = -128.e-6 0.0 -40.e-6 +slice.dom_hi = 128.e-6 0.0 0.96e-6 +slice.coarsening_ratio = 1 1 1 +slice.plot_int = -1 +slice.num_slice_snapshots_lab = 4 +slice.dt_slice_snapshots_lab = 3.3356409519815207e-12 +slice.particle_slice_width_lab = 2.e-6 diff --git a/Examples/Modules/charged_beam/inputs b/Examples/Modules/charged_beam/inputs deleted file mode 100644 index 18b645281..000000000 --- a/Examples/Modules/charged_beam/inputs +++ /dev/null @@ -1,50 +0,0 @@ -# Maximum number of time steps -max_step = 40 - -# number of grid points -amr.n_cell = 63 63 63 - -# Maximum allowable size of each subdomain in the problem domain; -# this is used to decompose the domain for parallel calculations. -amr.max_grid_size = 32 - -# 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 1 # Is periodic? -geometry.prob_lo = -20.e-6 -20.e-6 -20.e-6 # physical domain -geometry.prob_hi = 20.e-6 20.e-6 20.e-6 - -# Verbosity -warpx.verbose = 1 - -# Algorithms -algo.current_deposition = direct - -# CFL -warpx.cfl = 1.0 - -particles.nspecies = 1 -particles.species_names = electrons - -electrons.charge = -q_e -electrons.mass = m_e -electrons.injection_style = "NUniformPerCell" -electrons.num_particles_per_cell_each_dim = 2 2 2 - -electrons.xmin = -20.e-6 -electrons.xmax = 0.e-6 -electrons.ymin = -20.e-6 -electrons.ymax = 20.e-6 -electrons.zmin = -20.e-6 -electrons.zmax = 20.e-6 - -electrons.profile = constant -electrons.density = 1.e25 # number of electrons per m^3 - -electrons.momentum_distribution_type = "constant" -electrons.ux = 0.01 # ux = gamma*beta_x diff --git a/Examples/Modules/gaussian_beam/gaussian_beam_PICMI.py b/Examples/Modules/gaussian_beam/PICMI_inputs_gaussian_beam.py index a22e83794..ecc8f5a65 100644 --- a/Examples/Modules/gaussian_beam/gaussian_beam_PICMI.py +++ b/Examples/Modules/gaussian_beam/PICMI_inputs_gaussian_beam.py @@ -45,9 +45,9 @@ electrons = picmi.Species(particle_type='electron', name='electrons', initial_di protons = picmi.Species(particle_type='proton', name='protons', initial_distribution=proton_beam) sim = picmi.Simulation(solver = solver, - max_steps = 1000, + max_steps = 10, verbose = 1, - warpx_plot_int = 8, + warpx_plot_int = 10, warpx_current_deposition_algo = 'direct') sim.add_species(electrons, layout=picmi.PseudoRandomLayout(n_macroparticles=number_sim_particles)) diff --git a/Examples/Modules/gaussian_beam/inputs b/Examples/Modules/gaussian_beam/inputs deleted file mode 100644 index 46cd785f2..000000000 --- a/Examples/Modules/gaussian_beam/inputs +++ /dev/null @@ -1,100 +0,0 @@ -# Maximum number of time steps -max_step = 1000 - -# 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 = 8 # 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. -2. # physical domain -geometry.prob_hi = 2. 2. 2. - -# 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 = 2 -particles.species_names = electrons protons - -# -# The electron species information -# - -electrons.charge = -q_e -electrons.mass = m_e -electrons.injection_style = "gaussian_beam" -electrons.x_rms = 0.25 -electrons.y_rms = 0.25 -electrons.z_rms = 0.25 -electrons.x_m = 0. -electrons.y_m = 0. -electrons.z_m = 0. -electrons.npart = 32768 -electrons.q_tot = -8.010883097437485e-07 - -electrons.profile = "constant" -electrons.density = 1 -electrons.momentum_distribution_type = "radial_expansion" -electrons.u_over_r = -0.04 - -electrons.xmin = -2 -electrons.xmax = 2 -electrons.ymin = -2 -electrons.ymax = 2 -electrons.zmin = -2 -electrons.zmax = 2 - -# -# The proton species information -# - -protons.charge = q_e -protons.mass = m_p -protons.injection_style = "gaussian_beam" -protons.x_rms = 0.25 -protons.y_rms = 0.25 -protons.z_rms = 0.25 -protons.x_m = 0. -protons.y_m = 0. -protons.z_m = 0. -protons.npart = 32768 -protons.q_tot = 8.010883097437485e-07 - -protons.profile = "constant" -protons.density = 1 -protons.momentum_distribution_type = "radial_expansion" -protons.u_over_r = 0. - -protons.xmin = -2 -protons.xmax = 2 -protons.ymin = -2 -protons.ymax = 2 -protons.zmin = -2 -protons.zmax = 2 - -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 diff --git a/Examples/Modules/ionization/ionization_analysis.py b/Examples/Modules/ionization/analysis_ionization.py index f512eac6e..f512eac6e 100755 --- a/Examples/Modules/ionization/ionization_analysis.py +++ b/Examples/Modules/ionization/analysis_ionization.py diff --git a/Examples/Modules/laser_injection/Visualization.ipynb b/Examples/Modules/laser_injection/Visualization.ipynb deleted file mode 100644 index 68dd53ac7..000000000 --- a/Examples/Modules/laser_injection/Visualization.ipynb +++ /dev/null @@ -1,137 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Overview\n", - "\n", - "This a notebook that inspects the results of a WarpX simulation.\n", - "\n", - "# Instructions\n", - "\n", - "Execute the cells below one by one, by selecting them with your mouse and typing `Shift + Enter`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "# Import statements\n", - "import sys\n", - "from tqdm import tqdm\n", - "import yt, glob\n", - "yt.funcs.mylog.setLevel(50)\n", - "from IPython.display import clear_output\n", - "import numpy as np\n", - "from ipywidgets import interact, RadioButtons, IntSlider\n", - "import matplotlib.pyplot as plt\n", - "%matplotlib\n", - "\n", - "# Find iterations\n", - "file_list = glob.glob('plt?????')\n", - "iterations = [ int(file_name[3:]) for file_name in file_list ]" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Functions to plot the fields" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "def plot_field( iteration, field, slicing_direction='y', plotter='matplotlib' ):\n", - " ds = yt.load( './plt%05d/' %iteration )\n", - " all_data_level_0 = ds.covering_grid(level=0, \n", - " left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)\n", - " \n", - " if plotter == 'yt':\n", - " sl = yt.SlicePlot(ds, slicing_direction, field)\n", - " sl.set_log( field, False)\n", - " sl.annotate_grids()\n", - " # Show the new plot\n", - " clear_output()\n", - " sl.show()\n", - "\n", - " elif plotter == 'matplotlib':\n", - "\n", - " left_edge = ds.domain_left_edge.convert_to_mks()*1.e6\n", - " right_edge = ds.domain_right_edge.convert_to_mks()*1.e6\n", - " \n", - " if slicing_direction == 'x':\n", - " n = int( ds.domain_dimensions[0]//2 )\n", - " data2d = all_data_level_0[field][n, :, :]\n", - " extent = [ left_edge[2], right_edge[2], left_edge[1], right_edge[1] ]\n", - " elif slicing_direction == 'y':\n", - " n = int( ds.domain_dimensions[1]//2 )\n", - " data2d = all_data_level_0[field][:, n, :]\n", - " extent = [ left_edge[2], right_edge[2], left_edge[0], right_edge[0] ]\n", - " elif slicing_direction == 'z':\n", - " n = int( ds.domain_dimensions[2]//2 )\n", - " data2d = all_data_level_0[field][:, :, n]\n", - " extent = [ left_edge[1], right_edge[1], left_edge[0], right_edge[0] ]\n", - " plt.clf()\n", - " plt.title(\"%s at iteration %d\" %(field, iteration) )\n", - " plt.imshow( data2d, interpolation='nearest', cmap='viridis',\n", - " origin='lower', extent=extent, aspect='auto' )\n", - " plt.colorbar()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Interactive viewer" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "interact(plot_field, \n", - " iteration = IntSlider(min=min(iterations), max=max(iterations), step=iterations[1]-iterations[0]),\n", - " field = RadioButtons( options=['jx', 'jy', 'jz', 'Ex', 'Ey', 'Ez'], value='jz'),\n", - " slicing_direction = RadioButtons( options=[ 'x', 'y', 'z'], value='y'),\n", - " plotter = RadioButtons( options=['matplotlib', 'yt'] ) )" - ] - } - ], - "metadata": { - "anaconda-cloud": {}, - "kernelspec": { - "display_name": "Python [default]", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.5.2" - } - }, - "nbformat": 4, - "nbformat_minor": 1 -} diff --git a/Examples/Modules/laser_injection/laser_analysis.py b/Examples/Modules/laser_injection/analysis_laser.py index 1951bb29a..1951bb29a 100755 --- a/Examples/Modules/laser_injection/laser_analysis.py +++ b/Examples/Modules/laser_injection/analysis_laser.py diff --git a/Examples/Modules/nci_corrector/ncicorr_analysis.py b/Examples/Modules/nci_corrector/analysis_ncicorr.py index 439f40565..94dd2f838 100755 --- a/Examples/Modules/nci_corrector/ncicorr_analysis.py +++ b/Examples/Modules/nci_corrector/analysis_ncicorr.py @@ -11,11 +11,11 @@ fn = sys.argv[1] use_MR = re.search( 'nci_correctorMR', fn ) != None if use_MR: - energy_corrector_off = 1.e26 - energy_threshold = 3.e24 + energy_corrector_off = 5.e32 + energy_threshold = 1.e28 else: - energy_corrector_off = 1.e28 - energy_threshold = 5.e24 + energy_corrector_off = 1.5e26 + energy_threshold = 1.e24 # Check EB energy after 1000 timesteps filename = sys.argv[1] diff --git a/Examples/Modules/nci_corrector/inputs2d b/Examples/Modules/nci_corrector/inputs.2d index 0bcffbe85..61278a8e6 100644 --- a/Examples/Modules/nci_corrector/inputs2d +++ b/Examples/Modules/nci_corrector/inputs.2d @@ -14,8 +14,6 @@ amr.blocking_factor = 8 # Maximum level in hierarchy (for now must be 0, i.e., one level in total) amr.max_level = 0 -warpx.fine_tag_lo = -20.e-6 -20.e-6 -warpx.fine_tag_hi = 20.e-6 20.e-6 # Geometry geometry.coord_sys = 0 # 0: Cartesian diff --git a/Examples/Modules/qed/breit_wheeler/analysis_2d_tau_init.py b/Examples/Modules/qed/breit_wheeler/analysis_2d_tau_init.py new file mode 100755 index 000000000..850ecc0fe --- /dev/null +++ b/Examples/Modules/qed/breit_wheeler/analysis_2d_tau_init.py @@ -0,0 +1,31 @@ +#! /usr/bin/env python3 +import yt +import numpy as np +import scipy.stats as st +import sys + +# This script checks if photons initialized with Breit Wheeler process enabled +# do actually have an exponentially distributed optical depth + +# Tolerance +tol = 1e-2 + +def check(): + filename = sys.argv[1] + data_set = yt.load(filename) + + all_data = data_set.all_data() + res_tau = all_data["photons", 'particle_tau'] + + loc, scale = st.expon.fit(res_tau) + + # loc should be very close to 0, scale should be very close to 1 + assert(np.abs(loc - 0) < tol) + assert(np.abs(scale - 1) < tol) + +def main(): + check() + +if __name__ == "__main__": + main() + diff --git a/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init b/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init new file mode 100644 index 000000000..8e3dd29ca --- /dev/null +++ b/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init @@ -0,0 +1,80 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 1 +amr.n_cell = 128 128 +amr.max_grid_size = 128 # maximum size of each AMReX box, used to decompose the domain +amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain +amr.plot_int = 1 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 +amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported) + +################################# +############ NUMERICS ########### +################################# +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = energy-conserving +algo.particle_pusher = boris +interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z +interpolation.noy = 3 +interpolation.noz = 3 +warpx.verbose = 1 +warpx.do_dive_cleaning = 0 +warpx.plot_raw_fields = 0 +warpx.plot_raw_fields_guards = 0 +warpx.use_filter = 1 +warpx.cfl = 1. # if 1., the time step is set to its CFL limit +warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition +warpx.serialize_ics = 1 + +################################# +############ PLASMA ############# +################################# +particles.nspecies = 1 # number of species +particles.species_names = photons +particles.photon_species = photons +################################# + +photons.charge = -q_e +photons.mass = m_e +photons.injection_style = "NUniformPerCell" +photons.profile = "constant" +photons.xmin = -30e-6 +photons.ymin = -30e-6 +photons.zmin = -30e-6 +photons.xmax = 30e-6 +photons.ymax = 30e-6 +photons.zmax = 30e-6 +photons.num_particles_per_cell_each_dim = 2 2 +photons.density = 1e19 +photons.profile = "constant" +photons.momentum_distribution_type = "gaussian" +photons.ux_m = 0.0 +photons.uy_m = 0.0 +photons.uz_m = 0.0 +photons.ux_th = 100. +photons.uy_th = 100. +photons.uz_th = 100. +##########QED#################### +photons.do_qed = 1 +photons.do_qed_breit_wheeler = 1 +################################# + +#######QED ENGINE PARAMETERS##### +qed_bw.chi_min = 0.01 +qed_bw.ignore_tables_for_test = 1 +#qed_bw.generate_table = 1 +#qed_bw.tab_dndt_chi_min = 0.01 +#qed_bw.tab_dndt_chi_max = 100 +#qed_bw.tab_dndt_how_many = 200 +#qed_bw.tab_pair_chi_min = 0.01 +#qed_bw.tab_pair_chi_max = 100 +#qed_bw.tab_pair_chi_how_many = 30 +#qed_bw.tab_pair_frac_how_many = 30 +#qed_bw.save_table_in = "bw_table" +#qed_bw.load_table_from = "bw_table" +################################# diff --git a/Examples/Modules/qed/quantum_synchrotron/analysis_2d_tau_init.py b/Examples/Modules/qed/quantum_synchrotron/analysis_2d_tau_init.py new file mode 100755 index 000000000..05b313ee6 --- /dev/null +++ b/Examples/Modules/qed/quantum_synchrotron/analysis_2d_tau_init.py @@ -0,0 +1,36 @@ +#! /usr/bin/env python3 +import yt +import numpy as np +import scipy.stats as st +import sys + +# This script checks if electrons and positrons initialized with +# Quantum Synchrotron process enabled +# do actually have an exponentially distributed optical depth + +# Tolerance +tol = 1e-2 + +def check(): + filename = sys.argv[1] + data_set = yt.load(filename) + + all_data = data_set.all_data() + res_ele_tau = all_data["electrons", 'particle_tau'] + res_pos_tau = all_data["positrons", 'particle_tau'] + + loc_ele, scale_ele = st.expon.fit(res_ele_tau) + loc_pos, scale_pos = st.expon.fit(res_pos_tau) + + # loc should be very close to 0, scale should be very close to 1 + assert(np.abs(loc_ele - 0) < tol) + assert(np.abs(loc_pos - 0) < tol) + assert(np.abs(scale_ele - 1) < tol) + assert(np.abs(scale_pos - 1) < tol) + +def main(): + check() + +if __name__ == "__main__": + main() + diff --git a/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init b/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init new file mode 100644 index 000000000..03a2a5798 --- /dev/null +++ b/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init @@ -0,0 +1,105 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 1 +amr.n_cell = 128 128 +amr.max_grid_size = 128 # maximum size of each AMReX box, used to decompose the domain +amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain +amr.plot_int = 1 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 +amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported) + +################################# +############ NUMERICS ########### +################################# +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = energy-conserving +algo.particle_pusher = boris +interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z +interpolation.noy = 3 +interpolation.noz = 3 +warpx.verbose = 1 +warpx.do_dive_cleaning = 0 +warpx.plot_raw_fields = 0 +warpx.plot_raw_fields_guards = 0 +warpx.use_filter = 1 +warpx.cfl = 1. # if 1., the time step is set to its CFL limit +warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition +warpx.serialize_ics = 1 + + +################################# +############ PLASMA ############# +################################# +particles.nspecies = 2 # number of species +particles.species_names = electrons positrons +################################# + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.profile = "constant" +electrons.xmin = -30e-6 +electrons.ymin = -30e-6 +electrons.zmin = -30e-6 +electrons.xmax = 30e-6 +electrons.ymax = 30e-6 +electrons.zmax = 30e-6 +electrons.num_particles_per_cell_each_dim = 2 2 +electrons.density = 1e19 +electrons.profile = "constant" +electrons.momentum_distribution_type = "gaussian" +electrons.ux_m = 0.0 +electrons.uy_m = 0.0 +electrons.uz_m = 0.0 +electrons.ux_th = 100. +electrons.uy_th = 100. +electrons.uz_th = 100. +##########QED#################### +electrons.do_qed = 1 +electrons.do_qed_quantum_sync = 1 +################################# + +positrons.charge = q_e +positrons.mass = m_e +positrons.injection_style = "NUniformPerCell" +positrons.profile = "constant" +positrons.xmin = -30e-6 +positrons.ymin = -30e-6 +positrons.zmin = -30e-6 +positrons.xmax = 30e-6 +positrons.ymax = 30e-6 +positrons.zmax = 30e-6 +positrons.num_particles_per_cell_each_dim = 2 2 +positrons.density = 1e19 +positrons.profile = "constant" +positrons.momentum_distribution_type = "gaussian" +positrons.ux_m = 0.0 +positrons.uy_m = 0.0 +positrons.uz_m = 0.0 +positrons.ux_th = 100. +positrons.uy_th = 100. +positrons.uz_th = 100. +##########QED#################### +positrons.do_qed = 1 +positrons.do_qed_quantum_sync = 1 +################################# + +#######QED ENGINE PARAMETERS##### +qed_qs.chi_min = 0.001 +qed_qs.ignore_tables_for_test = 1 +#qed_qs.generate_table = 1 +#qed_qs.tab_dndt_chi_min = 0.001 +#qed_qs.tab_dndt_chi_max = 100 +#qed_qs.tab_dndt_how_many = 200 +#qed_qs.tab_em_chi_min = 0.01 +#qed_qs.tab_em_chi_max = 100 +#qed_qs.tab_em_chi_how_many = 2 +#qed_qs.tab_em_prob_how_many = 2 +#qed_qs.save_table_in = "qs_table" +#qed_qs.load_table_from = "qs_table" +################################# diff --git a/Examples/Modules/relativistic_space_charge_initialization/analysis.py b/Examples/Modules/relativistic_space_charge_initialization/analysis.py new file mode 100755 index 000000000..aad4534f5 --- /dev/null +++ b/Examples/Modules/relativistic_space_charge_initialization/analysis.py @@ -0,0 +1,75 @@ +#! /usr/bin/env python +""" +This script checks the space-charge initialization routine, by +verifying that the space-charge field of a Gaussian beam corresponds to +the expected theoretical field. +""" +import sys +import yt +import numpy as np +import matplotlib.pyplot as plt +import scipy.constants as scc +from scipy.special import gammainc +yt.funcs.mylog.setLevel(0) + +# Parameters from the Simulation +Qtot = -1.e-20 +r0 = 2.e-6 + +# Open data file +filename = sys.argv[1] +ds = yt.load( filename ) +# Extract data +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray().squeeze() +By_array = ad0['By'].to_ndarray() + +# Extract grid coordinates +Nx, Ny, Nz = ds.domain_dimensions +xmin, ymin, zmin = ds.domain_left_edge.v +Lx, Ly, Lz = ds.domain_width.v +x = xmin + Lx/Nx*(0.5+np.arange(Nx)) +y = ymin + Ly/Ny*(0.5+np.arange(Ny)) +z = zmin + Lz/Nz*(0.5+np.arange(Nz)) + +# Compute theoretical field +x_2d, y_2d, z_2d = np.meshgrid(x, y, z, indexing='ij') +r2 = x_2d**2 + y_2d**2 +factor = Qtot/scc.epsilon_0/(2*np.pi*r2) * (1-np.exp(-r2/(2*r0**2))) +factor_z = 1./(2*np.pi)**.5/r0 * np.exp(-z_2d**2/(2*r0**2)) +Ex_th = factor*factor_z*x_2d +Ey_th = factor*factor_z*y_2d + +# Plot theory and data +def make_2d(arr): + if arr.ndim == 3: + return arr[:,Ny//2,:] + else: + return arr +plt.figure(figsize=(10,10)) +plt.subplot(221) +plt.title('Ex: Theory') +plt.imshow(make_2d(Ex_th)) +plt.colorbar() +plt.subplot(222) +plt.title('Ex: Simulation') +plt.imshow(make_2d(Ex_array)) +plt.colorbar() +plt.subplot(223) +plt.title('By: Theory') +plt.imshow(make_2d(Ex_th/scc.c)) +plt.colorbar() +plt.subplot(224) +plt.title('Bz: Simulation') +plt.imshow(make_2d(By_array)) +plt.colorbar() + +plt.savefig('Comparison.png') + +# Automatically check the results +def check(E, E_th, label): + print( 'Relative error in %s: %.3f'%( + label, abs(E-E_th).max()/E_th.max())) + assert np.allclose( E, E_th, atol=0.1*E_th.max() ) + +check( Ex_array, Ex_th, 'Ex' ) diff --git a/Examples/Modules/relativistic_space_charge_initialization/inputs b/Examples/Modules/relativistic_space_charge_initialization/inputs new file mode 100644 index 000000000..eca38e074 --- /dev/null +++ b/Examples/Modules/relativistic_space_charge_initialization/inputs @@ -0,0 +1,30 @@ +max_step = 0 +amr.n_cell = 64 64 64 +amr.max_grid_size = 32 +amr.max_level = 0 + +amr.plot_int = 1 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 0 # Is periodic? +geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6 +geometry.prob_hi = 50.e-6 50.e-6 50.e-6 + +particles.nspecies = 1 +particles.species_names = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.initialize_self_fields = 1 +beam.self_fields_required_precision = 1.e-4 +beam.x_rms = 2.e-6 +beam.y_rms = 2.e-6 +beam.z_rms = 2.e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = 0.e-6 +beam.npart = 20000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 100.0 diff --git a/Examples/Modules/restart/inputs b/Examples/Modules/restart/inputs deleted file mode 100644 index 175c714b2..000000000 --- a/Examples/Modules/restart/inputs +++ /dev/null @@ -1,18 +0,0 @@ -# Basic simulation parameters -max_step = 500 -amr.n_cell = 256 256 -amr.max_grid_size = 32 -amr.max_level = 0 -particles.nspecies = 0 - -# Geometry parameters -geometry.coord_sys = 0 # 0: Cartesian -geometry.is_periodic = 0 0 0 # Is periodic? -geometry.prob_lo = -50.e-6 0. # physical domain -geometry.prob_hi = 50.e-6 30.e-6 - -# Number of iterations between consecutive checkpoint dumps -amr.check_int = 200 - -# Checkpoint file from which to restart the simulation. -# amr.restart = chk00400 diff --git a/Examples/Modules/space_charge_initialization/analysis.py b/Examples/Modules/space_charge_initialization/analysis.py new file mode 100755 index 000000000..67039bad1 --- /dev/null +++ b/Examples/Modules/space_charge_initialization/analysis.py @@ -0,0 +1,89 @@ +#! /usr/bin/env python +""" +This script checks the space-charge initialization routine, by +verifying that the space-charge field of a Gaussian beam corresponds to +the expected theoretical field. +""" +import sys +import yt +import numpy as np +import matplotlib.pyplot as plt +import scipy.constants as scc +from scipy.special import gammainc +yt.funcs.mylog.setLevel(0) + +# Parameters from the Simulation +Qtot = -1.e-20 +r0 = 2.e-6 + +# Open data file +filename = sys.argv[1] +ds = yt.load( filename ) +# Extract data +ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +Ex_array = ad0['Ex'].to_ndarray().squeeze() +if ds.dimensionality == 2: + # Rename the z dimension as y, so as to make this script work for 2d and 3d + Ey_array = ad0['Ez'].to_ndarray().squeeze() +elif ds.dimensionality == 3: + Ey_array = ad0['Ey'].to_ndarray() + Ez_array = ad0['Ez'].to_ndarray() + +# Extract grid coordinates +Nx, Ny, Nz = ds.domain_dimensions +xmin, ymin, zmin = ds.domain_left_edge.v +Lx, Ly, Lz = ds.domain_width.v +x = xmin + Lx/Nx*(0.5+np.arange(Nx)) +y = ymin + Ly/Ny*(0.5+np.arange(Ny)) +z = zmin + Lz/Nz*(0.5+np.arange(Nz)) + +# Compute theoretical field +if ds.dimensionality == 2: + x_2d, y_2d = np.meshgrid(x, y, indexing='ij') + r2 = x_2d**2 + y_2d**2 + factor = (Qtot/r0)/(2*np.pi*scc.epsilon_0*r2) * (1-np.exp(-r2/(2*r0**2))) + Ex_th = x_2d * factor + Ey_th = y_2d * factor +elif ds.dimensionality == 3: + x_2d, y_2d, z_2d = np.meshgrid(x, y, z, indexing='ij') + r2 = x_2d**2 + y_2d**2 + z_2d**2 + factor = Qtot/(4*np.pi*scc.epsilon_0*r2**1.5) * gammainc(3./2, r2/(2.*r0**2)) + Ex_th = factor*x_2d + Ey_th = factor*y_2d + Ez_th = factor*z_2d + +# Plot theory and data +def make_2d(arr): + if arr.ndim == 3: + return arr[:,:,Nz//2] + else: + return arr +plt.figure(figsize=(10,10)) +plt.subplot(221) +plt.title('Ex: Theory') +plt.imshow(make_2d(Ex_th)) +plt.colorbar() +plt.subplot(222) +plt.title('Ex: Simulation') +plt.imshow(make_2d(Ex_array)) +plt.colorbar() +plt.subplot(223) +plt.title('Ey: Theory') +plt.imshow(make_2d(Ey_th)) +plt.colorbar() +plt.subplot(224) +plt.title('Ey: Simulation') +plt.imshow(make_2d(Ey_array)) +plt.colorbar() +plt.savefig('Comparison.png') + +# Automatically check the results +def check(E, E_th, label): + print( 'Relative error in %s: %.3f'%( + label, abs(E-E_th).max()/E_th.max())) + assert np.allclose( E, E_th, atol=0.1*E_th.max() ) + +check( Ex_array, Ex_th, 'Ex' ) +check( Ey_array, Ey_th, 'Ey' ) +if ds.dimensionality == 3: + check( Ez_array, Ez_th, 'Ez' ) diff --git a/Examples/Modules/space_charge_initialization/inputs b/Examples/Modules/space_charge_initialization/inputs new file mode 100644 index 000000000..79309e585 --- /dev/null +++ b/Examples/Modules/space_charge_initialization/inputs @@ -0,0 +1,29 @@ +max_step = 0 +amr.n_cell = 64 64 64 +amr.max_grid_size = 32 +amr.max_level = 0 + +amr.plot_int = 1 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 0 # Is periodic? +geometry.prob_lo = -50.e-6 -50.e-6 -50.e-6 +geometry.prob_hi = 50.e-6 50.e-6 50.e-6 + +particles.nspecies = 1 +particles.species_names = beam +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.initialize_self_fields = 1 +beam.x_rms = 2.e-6 +beam.y_rms = 2.e-6 +beam.z_rms = 2.e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = 0.e-6 +beam.npart = 20000 +beam.q_tot = -1.e-20 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 0.0 |