diff options
| -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-- | Regression/WarpX-tests.ini | 33 | ||||
| -rw-r--r-- | Source/FieldSolver/WarpXPushFieldsEM.cpp | 64 | ||||
| -rw-r--r-- | Tools/read_raw_data.py | 2 |
8 files changed, 271 insertions, 129 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/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 0fcb88e05..b131249c8 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -102,6 +102,39 @@ analysisRoutine = Examples/Tests/PML/analysis_pml_ckc.py #doVis = 0 #analysisRoutine = Examples/Tests/PML/analysis_pml_psatd.py +[RigidInjection_lab] +buildDir = . +inputFile = Examples/Modules/RigidInjection/inputs.LabFrame +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Modules/RigidInjection/analysis_rigid_injection_LabFrame.py + +[RigidInjection_boost_backtransformed] +buildDir = . +inputFile = Examples/Modules/RigidInjection/inputs.BoostedFrame +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +outputFile = "lab_frame_data" +doComparison = 0 +aux1File = Tools/read_raw_data.py +analysisRoutine = Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py + [nci_corrector] buildDir = . inputFile = Examples/Modules/nci_corrector/inputs2d diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 1df05bc0f..e390dd61a 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -159,33 +159,29 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) auto const& Eyfab = Ey->array(mfi); auto const& Ezfab = Ez->array(mfi); if (do_nodal) { - amrex::ParallelFor(tbx, + amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bx_nodal(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); - }); - amrex::ParallelFor(tby, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_by_nodal(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bz_nodal(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy); }); } else if (WarpX::maxwell_fdtd_solver_id == 0) { - amrex::ParallelFor(tbx, + amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bx_yee(j,k,l,Bxfab,Eyfab,Ezfab,dtsdy,dtsdz); - }); - amrex::ParallelFor(tby, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_by_yee(j,k,l,Byfab,Exfab,Ezfab,dtsdx,dtsdz); - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bz_yee(j,k,l,Bzfab,Exfab,Eyfab,dtsdx,dtsdy,dxinv,xmin); @@ -198,23 +194,21 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - amrex::ParallelFor(tbx, + amrex::ParallelFor(tbx, tby, tbz, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bx_ckc(j,k,l,Bxfab,Eyfab,Ezfab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tby, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_by_ckc(j,k,l,Byfab,Exfab,Ezfab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_bz_ckc(j,k,l,Bzfab,Exfab,Eyfab, @@ -356,33 +350,29 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) auto const& jzfab = jz->array(mfi); if (do_nodal) { - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_nodal(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_nodal(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_nodal(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2); }); } else { - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_yee(j,k,l,Exfab,Byfab,Bzfab,jxfab,mu_c2_dt,dtsdy_c2,dtsdz_c2); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_yee(j,k,l,Eyfab,Bxfab,Bzfab,jyfab,mu_c2_dt,dtsdx_c2,dtsdz_c2,xmin); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_yee(j,k,l,Ezfab,Bxfab,Byfab,jzfab,mu_c2_dt,dtsdx_c2,dtsdy_c2,dxinv,xmin); @@ -393,17 +383,15 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { auto const& Ffab = F->array(mfi); if (WarpX::maxwell_fdtd_solver_id == 0) { - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_f_yee(j,k,l,Exfab,Ffab,dtsdx_c2); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_f_yee(j,k,l,Eyfab,Ffab,dtsdy_c2); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_f_yee(j,k,l,Ezfab,Ffab,dtsdz_c2); @@ -417,23 +405,21 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - amrex::ParallelFor(tex, + amrex::ParallelFor(tex, tey, tez, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ex_f_ckc(j,k,l,Exfab,Ffab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tey, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ey_f_ckc(j,k,l,Eyfab,Ffab, betaxy, betaxz, betayx, betayz, betazx, betazy, gammax, gammay, gammaz, alphax, alphay, alphaz); - }); - amrex::ParallelFor(tez, + }, [=] AMREX_GPU_DEVICE (int j, int k, int l) { warpx_push_ez_f_ckc(j,k,l,Ezfab,Ffab, @@ -616,7 +602,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // Rescale current in r-z mode since the inverse volume factor was not // included in the current deposition. - amrex::ParallelFor(tbr, + amrex::ParallelFor(tbr, tbt, tbz, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Wrap the current density deposited in the guard cells around @@ -630,8 +616,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // between on axis and off-axis factors const amrex::Real r = std::abs(rmin + (i - irmin + 0.5)*dr); Jr_arr(i,j,0) /= (2.*MathConst::pi*r); - }); - amrex::ParallelFor(tbt, + }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Wrap the current density deposited in the guard cells around @@ -649,8 +634,7 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu } else { Jt_arr(i,j,0) /= (2.*MathConst::pi*r); } - }); - amrex::ParallelFor(tbz, + }, [=] AMREX_GPU_DEVICE (int i, int j, int k) { // Wrap the current density deposited in the guard cells around diff --git a/Tools/read_raw_data.py b/Tools/read_raw_data.py index bca14b582..363c3a007 100644 --- a/Tools/read_raw_data.py +++ b/Tools/read_raw_data.py @@ -312,7 +312,7 @@ def _read_buffer(snapshot, header_fn, _component_names): for i in range(header.ncomp): comp_data = arr[i*size:(i+1)*size].reshape(shape, order='F') data = all_data[_component_names[i]] - data[[slice(l,h+1) for l, h in zip(lo, hi)]] = comp_data + data[tuple([slice(l,h+1) for l, h in zip(lo, hi)])] = comp_data all_data[_component_names[i]] = data return all_data |