From 4ef16e54c644d5018361ca836125237b052de954 Mon Sep 17 00:00:00 2001 From: Marco Garten Date: Thu, 22 Feb 2024 23:31:51 -0800 Subject: [PATCH] Improve Documentation for Laser-Ion Acceleration from Planar Target (#4569) * Update laser-ion acc example - add PICMI input - add vis script - update docs - add Python regression test - update phase space analysis script - add checksum regression test * Update docs with figures from hi-res runs * Add captions * Fix checksum regression tests * Change outputs to h5 for CI * Fix: Do Not Deselect Positions * Synchronize PICMI and Inputs More - Species selections - Filter details * Fix checksum regression test Update checksums in Python_LaserIonAcc2d Co-authored-by: Axel Huebl --- .../laser_ion/PICMI_inputs_2d.py | 313 ++++++++++++++++++ .../Physics_applications/laser_ion/README.rst | 71 +++- .../laser_ion/analysis_histogram_2D.py | 44 ++- .../Physics_applications/laser_ion/inputs_2d | 42 ++- .../Physics_applications/laser_ion/plot_2d.py | 229 +++++++++++++ .../benchmarks_json/LaserIonAcc2d.json | 51 +-- .../benchmarks_json/Python_LaserIonAcc2d.json | 34 ++ Regression/WarpX-tests.ini | 24 +- 8 files changed, 732 insertions(+), 76 deletions(-) create mode 100755 Examples/Physics_applications/laser_ion/PICMI_inputs_2d.py create mode 100644 Examples/Physics_applications/laser_ion/plot_2d.py create mode 100644 Regression/Checksum/benchmarks_json/Python_LaserIonAcc2d.json diff --git a/Examples/Physics_applications/laser_ion/PICMI_inputs_2d.py b/Examples/Physics_applications/laser_ion/PICMI_inputs_2d.py new file mode 100755 index 00000000000..844501992c3 --- /dev/null +++ b/Examples/Physics_applications/laser_ion/PICMI_inputs_2d.py @@ -0,0 +1,313 @@ +#!/usr/bin/env python3 + +from pywarpx import picmi + +# Physical constants +c = picmi.constants.c +q_e = picmi.constants.q_e + +# We only run 100 steps for tests +# Disable `max_step` below to run until the physical `stop_time`. +max_step = 100 +# time-scale with highly kinetic dynamics +stop_time = 0.2e-12 + +# proper resolution for 30 n_c (dx<=3.33nm) incl. acc. length +# (>=6x V100) +# --> choose larger `max_grid_size` and `blocking_factor` for 1 to 8 grids per GPU accordingly +#nx = 7488 +#nz = 14720 + +# Number of cells +nx = 384 +nz = 512 + +# Domain decomposition (deactivate `warpx_numprocs` in `picmi.Simulation` for this to take effect) +max_grid_size = 64 +blocking_factor = 32 + +# Physical domain +xmin = -7.5e-06 +xmax = 7.5e-06 +zmin = -5.0e-06 +zmax = 25.0e-06 + +# Create grid +grid = picmi.Cartesian2DGrid( + number_of_cells=[nx, nz], + lower_bound=[xmin, zmin], + upper_bound=[xmax, zmax], + lower_boundary_conditions=['open', 'open'], + upper_boundary_conditions=['open', 'open'], + lower_boundary_conditions_particles=['absorbing', 'absorbing'], + upper_boundary_conditions_particles=['absorbing', 'absorbing'], + warpx_max_grid_size=max_grid_size, + warpx_blocking_factor=blocking_factor) + +# Particles: plasma parameters +# critical plasma density +nc = 1.742e27 # [m^-3] 1.11485e21 * 1.e6 / 0.8**2 +# number density: "fully ionized" electron density as reference +# [material 1] cryogenic H2 +n0 = 30.0 # [n_c] +# [material 2] liquid crystal +# n0 = 192 +# [material 3] PMMA +# n0 = 230 +# [material 4] Copper (ion density: 8.49e28/m^3; times ionization level) +# n0 = 1400 +plasma_density = n0 * nc +preplasma_L = 0.05e-6 # [m] scale length (>0) +preplasma_Lcut = 2.0e-6 # [m] hard cutoff from surface +plasma_r0 = 2.5e-6 # [m] radius or half-thickness +plasma_eps_z = 0.05e-6 # [m] small offset in z to make zmin, zmax interval larger than 2*(r0 + Lcut) +plasma_creation_limit_z = plasma_r0 + preplasma_Lcut + plasma_eps_z # [m] upper limit in z for particle creation + +plasma_xmin = None +plasma_ymin = None +plasma_zmin = -plasma_creation_limit_z +plasma_xmax = None +plasma_ymax = None +plasma_zmax = plasma_creation_limit_z + +density_expression_str = f'{plasma_density}*((abs(z)<={plasma_r0}) + (abs(z)<{plasma_r0}+{preplasma_Lcut}) * (abs(z)>{plasma_r0}) * exp(-(abs(z)-{plasma_r0})/{preplasma_L}))' + +slab_with_ramp_dist_hydrogen = picmi.AnalyticDistribution( + density_expression=density_expression_str, + lower_bound=[plasma_xmin, plasma_ymin, plasma_zmin], + upper_bound=[plasma_xmax, plasma_ymax, plasma_zmax] +) + +# thermal velocity spread for electrons in gamma*beta +ux_th = .01 +uz_th = .01 + +slab_with_ramp_dist_electrons = picmi.AnalyticDistribution( + density_expression=density_expression_str, + lower_bound=[plasma_xmin, plasma_ymin, plasma_zmin], + upper_bound=[plasma_xmax, plasma_ymax, plasma_zmax], + # if `momentum_expressions` and `momentum_spread_expressions` are unset, + # a Gaussian momentum distribution is assumed given that `rms_velocity` has any non-zero elements + rms_velocity=[c*ux_th, 0., c*uz_th] # thermal velocity spread in m/s +) + +# TODO: add additional attributes orig_x and orig_z +electrons = picmi.Species( + particle_type='electron', + name='electrons', + initial_distribution=slab_with_ramp_dist_electrons, +) + +# TODO: add additional attributes orig_x and orig_z +hydrogen = picmi.Species( + particle_type='proton', + name='hydrogen', + initial_distribution=slab_with_ramp_dist_hydrogen +) + +# Laser +# e_max = a0 * 3.211e12 / lambda_0[mu] +# a0 = 16, lambda_0 = 0.8mu -> e_max = 64.22 TV/m +e_max = 64.22e12 +position_z = -4.0e-06 +profile_t_peak = 50.e-15 +profile_focal_distance = 4.0e-06 +laser = picmi.GaussianLaser( + wavelength=0.8e-06, + waist=4.e-06, + duration=30.e-15, + focal_position=[0, 0, profile_focal_distance + position_z], + centroid_position=[0, 0, position_z - c * profile_t_peak], + propagation_direction=[0, 0, 1], + polarization_direction=[1, 0, 0], + E0=e_max, + fill_in=False) +laser_antenna = picmi.LaserAntenna( + position=[0., 0., position_z], + normal_vector=[0, 0, 1]) + +# Electromagnetic solver +solver = picmi.ElectromagneticSolver( + grid=grid, + method='Yee', + cfl=0.999, + divE_cleaning=0, + #warpx_pml_ncell=10 +) + +# Diagnostics +particle_diag = picmi.ParticleDiagnostic( + name='Python_LaserIonAcc2d_plt', + period=100, + write_dir='./diags', + warpx_format='openpmd', + warpx_openpmd_backend='h5', + # demonstration of a spatial and momentum filter + warpx_plot_filter_function='(uz>=0) * (x<1.0e-6) * (x>-1.0e-6)' +) +# reduce resolution of output fields +coarsening_ratio = [4, 4] +ncell_field = [] +for (ncell_comp, cr) in zip([nx,nz], coarsening_ratio): + ncell_field.append(int(ncell_comp/cr)) +field_diag = picmi.FieldDiagnostic( + name='Python_LaserIonAcc2d_plt', + grid=grid, + period=100, + number_of_cells=ncell_field, + data_list=['B', 'E', 'J', 'rho', 'rho_electrons', 'rho_hydrogen'], + write_dir='./diags', + warpx_format='openpmd', + warpx_openpmd_backend='h5' +) + +particle_fw_diag = picmi.ParticleDiagnostic( + name='openPMDfw', + period=100, + write_dir='./diags', + warpx_format='openpmd', + warpx_openpmd_backend='h5', + warpx_plot_filter_function='(uz>=0) * (x<1.0e-6) * (x>-1.0e-6)' +) + +particle_bw_diag = picmi.ParticleDiagnostic( + name='openPMDbw', + period=100, + write_dir='./diags', + warpx_format='openpmd', + warpx_openpmd_backend='h5', + warpx_plot_filter_function='(uz<0)' +) + +# histograms with 2.0 degree acceptance angle in fw direction +# 2 deg * pi / 180 : 0.03490658503 rad +# half-angle +/- : 0.017453292515 rad +histuH_rdiag = picmi.ReducedDiagnostic( + diag_type='ParticleHistogram', + name='histuH', + period=100, + species=hydrogen, + bin_number=1000, + bin_min=0.0, + bin_max=0.474, # 100 MeV protons + histogram_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)', + filter_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)') + +histue_rdiag = picmi.ReducedDiagnostic( + diag_type='ParticleHistogram', + name='histue', + period=100, + species=electrons, + bin_number=1000, + bin_min=0.0, + bin_max=197.0, # 100 MeV electrons + histogram_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)', + filter_function='u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)') + +# just a test entry to make sure that the histogram filter is purely optional: +# this one just records uz of all hydrogen ions, independent of their pointing +histuzAll_rdiag = picmi.ReducedDiagnostic( + diag_type='ParticleHistogram', + name='histuzAll', + period=100, + species=hydrogen, + bin_number=1000, + bin_min=-0.474, + bin_max=0.474, + histogram_function='uz') + +field_probe_z_rdiag = picmi.ReducedDiagnostic( + diag_type='FieldProbe', + name='FieldProbe_Z', + period=100, + integrate=0, + probe_geometry='Line', + x_probe=0.0, + z_probe=-5.0e-6, + x1_probe=0.0, + z1_probe=25.0e-6, + resolution=3712) + +field_probe_scat_point_rdiag = picmi.ReducedDiagnostic( + diag_type='FieldProbe', + name='FieldProbe_ScatPoint', + period=1, + integrate=0, + probe_geometry='Point', + x_probe=0.0, + z_probe=15.0e-6) + +field_probe_scat_line_rdiag = picmi.ReducedDiagnostic( + diag_type='FieldProbe', + name='FieldProbe_ScatLine', + period=100, + integrate=1, + probe_geometry='Line', + x_probe=-2.5e-6, + z_probe=15.0e-6, + x1_probe=2.5e-6, + z1_probe=15e-6, + resolution=201) + +load_balance_costs_rdiag = picmi.ReducedDiagnostic( + diag_type='LoadBalanceCosts', + name='LBC', + period=100) + +# Set up simulation +sim = picmi.Simulation( + solver=solver, + max_time=stop_time, # need to remove `max_step` to run this far + verbose=1, + particle_shape='cubic', + warpx_numprocs=[1, 2], # deactivate `numprocs` for dynamic load balancing + warpx_use_filter=1, + warpx_load_balance_intervals=100, + warpx_load_balance_costs_update='heuristic' +) + +# Add plasma electrons +sim.add_species( + electrons, + layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[2,2]) + # for more realistic simulations, try to avoid that macro-particles represent more than 1 n_c + #layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[4,8]) +) + +# Add hydrogen ions +sim.add_species( + hydrogen, + layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[2,2]) + # for more realistic simulations, try to avoid that macro-particles represent more than 1 n_c + #layout=picmi.GriddedLayout(grid=grid, n_macroparticle_per_cell=[4,8]) +) + +# Add laser +sim.add_laser( + laser, + injection_method=laser_antenna) + +# Add full diagnostics +sim.add_diagnostic(particle_diag) +sim.add_diagnostic(field_diag) +sim.add_diagnostic(particle_fw_diag) +sim.add_diagnostic(particle_bw_diag) +# Add reduced diagnostics +sim.add_diagnostic(histuH_rdiag) +sim.add_diagnostic(histue_rdiag) +sim.add_diagnostic(histuzAll_rdiag) +sim.add_diagnostic(field_probe_z_rdiag) +sim.add_diagnostic(field_probe_scat_point_rdiag) +sim.add_diagnostic(field_probe_scat_line_rdiag) +sim.add_diagnostic(load_balance_costs_rdiag) +# TODO: make ParticleHistogram2D available + +# Write input file that can be used to run with the compiled version +sim.write_input_file(file_name='inputs_2d_picmi') + +# Initialize inputs and WarpX instance +sim.initialize_inputs() +sim.initialize_warpx() + +# Advance simulation until last time step +sim.step(max_step) diff --git a/Examples/Physics_applications/laser_ion/README.rst b/Examples/Physics_applications/laser_ion/README.rst index f382b36bc85..18a427aa1c7 100644 --- a/Examples/Physics_applications/laser_ion/README.rst +++ b/Examples/Physics_applications/laser_ion/README.rst @@ -4,18 +4,15 @@ Laser-Ion Acceleration with a Planar Target =========================================== This example shows how to model laser-ion acceleration with planar targets of solid density :cite:p:`ex-Wilks2001,ex-Bulanov2008,ex-Macchi2013`. -The acceleration mechanism in this scenario depends on target parameters +The acceleration mechanism in this scenario depends on target parameters. Although laser-ion acceleration requires full 3D modeling for adequate description of the acceleration dynamics, especially the acceleration field lengths and decay times, this example models a 2D example. 2D modeling can often hint at a qualitative overview of the dynamics, but mostly saves computational costs since the plasma frequency (and Debye length) of the plasma determines the resolution need in laser-solid interaction modeling. -.. note:: - - TODO: The Python (PICMI) input file needs to be created. - .. note:: The resolution of this 2D case is extremely low by default. + This includes spatial and temporal resolution, but also the number of macro-particles per cell representing the target density for proper phase space sampling. You will need a computing cluster for adequate resolution of the target density, see comments in the input file. .. warning:: @@ -30,18 +27,19 @@ Run This example can be run **either** as: -* **Python** script: (*TODO*) or -* WarpX **executable** using an input file: ``warpx.2d inputs_2d`` +* **Python** script: ``mpiexec -n 2 python3 PICMI_inputs_2d.py`` or +* WarpX **executable** using an input file: ``mpiexec -n 2 warpx.2d inputs_2d`` -For `MPI-parallel `__ runs, prefix these lines with ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system. +For `MPI-parallel `__ runs on computing clusters, change the prefix to ``mpiexec -n ...`` or ``srun -n ...``, depending on the system and number of MPI ranks you want to allocate. .. tab-set:: .. tab-item:: Python: Script - .. note:: + .. literalinclude:: PICMI_inputs_2d.py + :language: python3 + :caption: You can copy this file from ``Examples/Physics_applications/laser_ion/PICMI_inputs_2d.py``. - TODO: This input file should be created following the ``inputs_2d`` file. .. tab-item:: Executable: Input File @@ -52,14 +50,61 @@ For `MPI-parallel `__ runs, prefix these lines with ` Analyze ------- -.. note:: +.. _fig-tnsa-ps-electrons-pinhole: + +.. figure:: https://user-images.githubusercontent.com/5416860/295003882-c755fd47-4bb3-4439-9319-c48214cbaafd.png + :alt: Longitudinal phase space of forward-flying electrons in a 2 degree opening angle. + :width: 100% + + Longitudinal phase space of forward-flying electrons in a 2 degree opening angle. + +.. _fig-tnsa-ps-protons-pinhole: + +.. figure:: https://user-images.githubusercontent.com/5416860/295003988-dea3dfb7-0d55-4616-b32d-061fb429f9ac.png + :alt: Longitudinal phase space of forward-flying protons in a 2 degree opening angle. + :width: 100% + + Longitudinal phase space of forward-flying protons in a 2 degree opening angle. + +Time-resolved phase electron space analysis as in :numref:`fig-tnsa-ps-electrons-pinhole` gives information about, e.g., how laser energy is locally converted into electron kinetic energy. +Later in time, ion phase spaces like :numref:`fig-tnsa-ps-protons-pinhole` can reveal where accelerated ion populations originate. - This section is TODO. +.. dropdown:: Script ``analysis_histogram_2D.py`` + .. literalinclude:: analysis_histogram_2D.py + :language: python3 + :caption: You can copy this file from ``Examples/Physics_applications/laser_ion/analysis_histogram_2D.py``. Visualize --------- .. note:: - This section is TODO. + The following images for densities and electromagnetic fields were created with a run on 64 NVidia A100 GPUs featuring a total number of cells of ``nx = 8192`` and ``nz = 16384``, as well as 64 particles per cell per species. + +.. _fig-tnsa-densities: + +.. figure:: https://user-images.githubusercontent.com/5416860/296338802-8059c39c-0be8-4e4d-b41b-f976b626bd7f.png + :alt: Particle densities for electrons (top), protons (middle), and electrons again in logarithmic scale (bottom). + :width: 80% + + Particle densities for electrons (top), protons (middle), and electrons again in logarithmic scale (bottom). + +Particle density output illustrates the evolution of the target in time and space. +Logarithmic scales can help to identify where the target becomes transparent for the laser pulse (bottom panel in :numref:`fig-tnsa-densities` ). + +.. _fig-tnsa-fields: + +.. figure:: https://user-images.githubusercontent.com/5416860/296338609-a49eee7f-6793-4b55-92f1-0b887e6437ab.png + :alt: Electromagnetic field visualization for E_x (top), B_y (middle), and E_z (bottom). + :width: 80% + + Electromagnetic field visualization for :math:`E_x` (top), :math:`B_y` (middle), and :math:`E_z` (bottom). + +Electromagnetic field output shows where the laser field is strongest at a given point in time, and where accelerating fields build up :numref:`fig-tnsa-fields`. + +.. dropdown:: Script ``plot_2d.py`` + + .. literalinclude:: plot_2d.py + :language: python3 + :caption: You can copy this file from ``Examples/Physics_applications/laser_ion/plot_2d.py``. diff --git a/Examples/Physics_applications/laser_ion/analysis_histogram_2D.py b/Examples/Physics_applications/laser_ion/analysis_histogram_2D.py index 175d8eed568..a262a2373e5 100644 --- a/Examples/Physics_applications/laser_ion/analysis_histogram_2D.py +++ b/Examples/Physics_applications/laser_ion/analysis_histogram_2D.py @@ -6,6 +6,7 @@ import matplotlib.colors as colors import matplotlib.pyplot as plt +import numpy as np from openpmd_viewer import OpenPMDTimeSeries parser = argparse.ArgumentParser(description='Process a 2D histogram name and an integer.') @@ -18,30 +19,43 @@ ts = OpenPMDTimeSeries(path) it = ts.iterations -data, info = ts.get_field(field = "data", iteration = 0, plot = True) +data, info = ts.get_field(field="data", iteration=0, plot=True) print('The available iterations of the simulation are:', it) -print('The axis of the histogram are (0: ordinate ; 1: abscissa):', info.axes) +print('The axes of the histogram are (0: ordinate ; 1: abscissa):', info.axes) print('The data shape is:', data.shape) # Add the simulation time to the title once this information # is available in the "info" FieldMetaInformation object. if args.iter == 'All' : - for i in it : + for it_idx, i in enumerate(it): plt.figure() - data, info = ts.get_field(field = "data", iteration = i, plot = False) - plt.imshow(data, aspect="auto", norm=colors.LogNorm(), extent=info.imshow_extent) - plt.title(args.hist2D + " (iteration %d)" % i) - plt.xlabel(info.axes[1] + ' (m)') - plt.ylabel(info.axes[0] + ' (m.s-1)') + data, info = ts.get_field(field="data", iteration=i, plot=False) + abscissa_name = info.axes[1] # This might be 'z' or something else + abscissa_values = getattr(info, abscissa_name, None) + ordinate_name = info.axes[0] # This might be 'z' or something else + ordinate_values = getattr(info, ordinate_name, None) + + plt.pcolormesh(abscissa_values/1e-6, ordinate_values, data, norm=colors.LogNorm(), rasterized=True) + plt.title(args.hist2D + f" Time: {ts.t[it_idx]:.2e} s (Iteration: {i:d})") + plt.xlabel(info.axes[1]+r' ($\mu$m)') + plt.ylabel(info.axes[0]+r' ($m_\mathrm{species} c$)') plt.colorbar() - plt.savefig('Histogram_2D_' + args.hist2D + '_iteration_' + str(i) + '.pdf') + plt.tight_layout() + plt.savefig('Histogram_2D_' + args.hist2D + '_iteration_' + str(i) + '.png') else : i = int(args.iter) + it_idx = np.where(i == it)[0][0] plt.figure() - data, info = ts.get_field(field = "data", iteration = i, plot = False) - plt.imshow(data, aspect="auto", norm=colors.LogNorm(), extent=info.imshow_extent) - plt.title(args.hist2D + " (iteration %d)" % i) - plt.xlabel(info.axes[1] + ' (m)') - plt.ylabel(info.axes[0] + ' (m.s-1)') + data, info = ts.get_field(field="data", iteration=i, plot=False) + abscissa_name = info.axes[1] # This might be 'z' or something else + abscissa_values = getattr(info, abscissa_name, None) + ordinate_name = info.axes[0] # This might be 'z' or something else + ordinate_values = getattr(info, ordinate_name, None) + + plt.pcolormesh(abscissa_values/1e-6, ordinate_values, data, norm=colors.LogNorm(), rasterized=True) + plt.title(args.hist2D + f" Time: {ts.t[it_idx]:.2e} s (Iteration: {i:d})") + plt.xlabel(info.axes[1]+r' ($\mu$m)') + plt.ylabel(info.axes[0]+r' ($m_\mathrm{species} c$)') plt.colorbar() - plt.savefig('Histogram_2D_' + args.hist2D + '_iteration_' + str(i) + '.pdf') + plt.tight_layout() + plt.savefig('Histogram_2D_' + args.hist2D + '_iteration_' + str(i) + '.png') diff --git a/Examples/Physics_applications/laser_ion/inputs_2d b/Examples/Physics_applications/laser_ion/inputs_2d index 8fe5acf5f46..5ad8334e9ef 100644 --- a/Examples/Physics_applications/laser_ion/inputs_2d +++ b/Examples/Physics_applications/laser_ion/inputs_2d @@ -2,6 +2,9 @@ # Domain, Resolution & Numerics # +# We only run 100 steps for tests +# Disable `max_step` below to run until the physical `stop_time`. +max_step = 100 # time-scale with highly kinetic dynamics stop_time = 0.2e-12 # [s] # time-scale for converged ion energy @@ -9,12 +12,12 @@ stop_time = 0.2e-12 # [s] # - ions will start to leave the box #stop_time = 1.0e-12 # [s] -# quick tests at ultra-low res. (CI) -#amr.n_cell = 384 512 +# quick tests at ultra-low res. (for CI, and local computer) +amr.n_cell = 384 512 # proper resolution for 10 n_c excl. acc. length # (>=1x V100) -amr.n_cell = 2688 3712 +#amr.n_cell = 2688 3712 # proper resolution for 30 n_c (dx<=3.33nm) incl. acc. length # (>=6x V100) @@ -94,6 +97,8 @@ particles.species_names = electrons hydrogen hydrogen.species_type = hydrogen hydrogen.injection_style = NUniformPerCell hydrogen.num_particles_per_cell_each_dim = 2 2 +# for more realistic simulations, try to avoid that macro-particles represent more than 1 n_c +#hydrogen.num_particles_per_cell_each_dim = 4 8 hydrogen.momentum_distribution_type = at_rest # minimum and maximum z position between which particles are initialized # --> should be set for dense targets limit memory consumption during initialization @@ -107,6 +112,8 @@ hydrogen.attribute.orig_z(x,y,z,ux,uy,uz,t) = "z" electrons.species_type = electron electrons.injection_style = NUniformPerCell electrons.num_particles_per_cell_each_dim = 2 2 +# for more realistic simulations, try to avoid that macro-particles represent more than 1 n_c +#electrons.num_particles_per_cell_each_dim = 4 8 electrons.momentum_distribution_type = "gaussian" electrons.ux_th = .01 electrons.uz_th = .01 @@ -200,22 +207,19 @@ diag1.diag_type = Full diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen # reduce resolution of output fields diag1.coarsening_ratio = 4 4 -diag1.electrons.variables = w ux uy uz -diag1.hydrogen.variables = w ux uz orig_x orig_z # demonstration of a spatial and momentum filter diag1.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6) diag1.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6) +diag1.format = openpmd +diag1.openpmd_backend = h5 openPMDfw.intervals = 100 openPMDfw.diag_type = Full openPMDfw.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen # reduce resolution of output fields openPMDfw.coarsening_ratio = 4 4 -openPMDfw.electrons.variables = w ux uy uz -openPMDfw.hydrogen.variables = w ux uy uz orig_x orig_z openPMDfw.format = openpmd openPMDfw.openpmd_backend = h5 -openPMDfw.species = electrons hydrogen # demonstration of a spatial and momentum filter openPMDfw.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6) openPMDfw.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6) @@ -225,11 +229,8 @@ openPMDbw.diag_type = Full openPMDbw.fields_to_plot = rho_hydrogen # reduce resolution of output fields openPMDbw.coarsening_ratio = 4 4 -openPMDbw.electrons.variables = w ux uy uz -openPMDbw.hydrogen.variables = w ux uy uz orig_x orig_z openPMDbw.format = openpmd openPMDbw.openpmd_backend = h5 -openPMDbw.species = electrons hydrogen # demonstration of a momentum filter openPMDbw.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz<0) openPMDbw.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz<0) @@ -246,21 +247,19 @@ warpx.reduced_diags_names = histuH histue histuzAll FieldProbe histuH.type = ParticleHistogram histuH.intervals = 100 -histuH.path = "./" histuH.species = hydrogen histuH.bin_number = 1000 histuH.bin_min = 0.0 -histuH.bin_max = 35.0 +histuH.bin_max = 0.474 # 100 MeV protons histuH.histogram_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)" histuH.filter_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)" histue.type = ParticleHistogram histue.intervals = 100 -histue.path = "./" histue.species = electrons histue.bin_number = 1000 histue.bin_min = 0.0 -histue.bin_max = 0.1 +histue.bin_max = 197 # 100 MeV electrons histue.histogram_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, sqrt(u2), 0.0)" histue.filter_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(acos(uz / sqrt(u2))) <= 0.017453, 0)" @@ -268,11 +267,10 @@ histue.filter_function(t,x,y,z,ux,uy,uz) = "u2=ux*ux+uy*uy+uz*uz; if(u2>0, abs(a # this one just records uz of all hydrogen ions, independent of their pointing histuzAll.type = ParticleHistogram histuzAll.intervals = 100 -histuzAll.path = "./" histuzAll.species = hydrogen histuzAll.bin_number = 1000 -histuzAll.bin_min = -35.0 -histuzAll.bin_max = 35.0 +histuzAll.bin_min = -0.474 +histuzAll.bin_max = 0.474 histuzAll.histogram_function(t,x,y,z,ux,uy,uz) = "uz" FieldProbe_Z.type = FieldProbe @@ -313,8 +311,8 @@ PhaseSpaceIons.bin_number_abs = 1000 PhaseSpaceIons.bin_number_ord = 1000 PhaseSpaceIons.bin_min_abs = -5.e-6 PhaseSpaceIons.bin_max_abs = 25.e-6 -PhaseSpaceIons.bin_min_ord = -20.e-6 -PhaseSpaceIons.bin_max_ord = 20.e-6 +PhaseSpaceIons.bin_min_ord = -0.474 +PhaseSpaceIons.bin_max_ord = 0.474 PhaseSpaceIons.histogram_function_abs(t,x,y,z,ux,uy,uz,w) = "z" PhaseSpaceIons.histogram_function_ord(t,x,y,z,ux,uy,uz,w) = "uz" PhaseSpaceIons.value_function(t,x,y,z,ux,uy,uz,w) = "w" @@ -327,8 +325,8 @@ PhaseSpaceElectrons.bin_number_abs = 1000 PhaseSpaceElectrons.bin_number_ord = 1000 PhaseSpaceElectrons.bin_min_abs = -5.e-6 PhaseSpaceElectrons.bin_max_abs = 25.e-6 -PhaseSpaceElectrons.bin_min_ord = -20 -PhaseSpaceElectrons.bin_max_ord = 20 +PhaseSpaceElectrons.bin_min_ord = -197 +PhaseSpaceElectrons.bin_max_ord = 197 PhaseSpaceElectrons.histogram_function_abs(t,x,y,z,ux,uy,uz,w) = "z" PhaseSpaceElectrons.histogram_function_ord(t,x,y,z,ux,uy,uz,w) = "uz" PhaseSpaceElectrons.value_function(t,x,y,z,ux,uy,uz,w) = "w" diff --git a/Examples/Physics_applications/laser_ion/plot_2d.py b/Examples/Physics_applications/laser_ion/plot_2d.py new file mode 100644 index 00000000000..e85782a7d23 --- /dev/null +++ b/Examples/Physics_applications/laser_ion/plot_2d.py @@ -0,0 +1,229 @@ +#!/usr/bin/env python3 + +# Copyright 2023 The WarpX Community +# +# This file is part of WarpX. +# +# Authors: Marco Garten +# License: BSD-3-Clause-LBNL +# +# This script plots the densities and fields of a 2D laser-ion acceleration simulation. + + +import argparse +import os +import re + +from matplotlib.colors import TwoSlopeNorm +import matplotlib.pyplot as plt +import numpy as np +from openpmd_viewer import OpenPMDTimeSeries +import pandas as pd +import scipy.constants as sc + +plt.rcParams.update({'font.size':16}) + +def create_analysis_dir(directory): + if not os.path.exists(directory): + os.makedirs(directory) + + +def visualize_density_iteration(ts, iteration, out_dir): + """ + Visualize densities and fields of a single iteration. + + :param ts: OpenPMDTimeSeries + :param iteration: Output iteration (simulation timestep) + :param out_dir: Directory for PNG output + :return: + """ + # Physics parameters + lambda_L = 800e-9 # Laser wavelength in meters + omega_L = 2 * np.pi * sc.c / lambda_L # Laser frequency in seconds + n_c = sc.m_e * sc.epsilon_0 * omega_L**2 / sc.elementary_charge**2 # Critical plasma density in meters^(-3) + micron = 1e-6 + + # Simulation parameters + n_e0 = 30 + n_max = 2 * n_e0 + nr = 1 # Number to decrease resolution + + # Data fetching + it = iteration + ii = np.where(ts.iterations == it)[0][0] + + time = ts.t[ii] + rho_e, rho_e_info = ts.get_field(field="rho_electrons", iteration=it) + rho_d, rho_d_info = ts.get_field(field="rho_hydrogen", iteration=it) + + # Rescale to critical density + rho_e = rho_e / (sc.elementary_charge * n_c) + rho_d = rho_d / (sc.elementary_charge * n_c) + + # Axes setup + fig, axs = plt.subplots(3, 1, figsize=(5, 8)) + xax, zax = rho_e_info.x, rho_e_info.z + + # Plotting + # Electron density + im0 = axs[0].pcolormesh(zax[::nr]/micron, xax[::nr]/micron, -rho_e.T[::nr, ::nr], + vmin=0, vmax=n_max, cmap="Reds", rasterized=True) + plt.colorbar(im0, ax=axs[0], label=r"$n_\mathrm{\,e}\ (n_\mathrm{c})$") + + # Hydrogen density + im1 = axs[1].pcolormesh(zax[::nr]/micron, xax[::nr]/micron, rho_d.T[::nr, ::nr], + vmin=0, vmax=n_max, cmap="Blues", rasterized=True) + plt.colorbar(im1, ax=axs[1], label=r"$n_\mathrm{\,H}\ (n_\mathrm{c})$") + + # Masked electron density + divnorm = TwoSlopeNorm(vmin=-7., vcenter=0., vmax=2) + masked_data = np.ma.masked_where(rho_e.T == 0, rho_e.T) + my_cmap = plt.cm.PiYG_r.copy() + my_cmap.set_bad(color='black') + im2 = axs[2].pcolormesh(zax[::nr]/micron, xax[::nr]/micron, np.log(-masked_data[::nr, ::nr]), + norm=divnorm, cmap=my_cmap, rasterized=True) + plt.colorbar(im2, ax=axs[2], ticks=[-6, -3, 0, 1, 2], extend='both', + label=r"$\log n_\mathrm{\,e}\ (n_\mathrm{c})$") + + # Axis labels and title + for ax in axs: + ax.set_aspect(1.0) + ax.set_ylabel(r"$x$ ($\mu$m)") + for ax in axs[:-1]: + ax.set_xticklabels([]) + axs[2].set_xlabel(r"$z$ ($\mu$m)") + fig.suptitle(f"Iteration: {it}, Time: {time/1e-15:.1f} fs") + + plt.tight_layout() + + plt.savefig(f"{out_dir}/densities_{it:06d}.png") + +def visualize_field_iteration(ts, iteration, out_dir): + + # Additional parameters + nr = 1 # Number to decrease resolution + micron = 1e-6 + + # Data fetching + it = iteration + ii = np.where(ts.iterations == it)[0][0] + time = ts.t[ii] + + Ex, Ex_info = ts.get_field(field="E", coord="x", iteration=it) + Exmax = np.max(np.abs([np.min(Ex),np.max(Ex)])) + By, By_info = ts.get_field(field="B", coord="y", iteration=it) + Bymax = np.max(np.abs([np.min(By),np.max(By)])) + Ez, Ez_info = ts.get_field(field="E", coord="z", iteration=it) + Ezmax = np.max(np.abs([np.min(Ez),np.max(Ez)])) + + # Axes setup + fig,axs = plt.subplots(3, 1, figsize=(5, 8)) + xax, zax = Ex_info.x, Ex_info.z + + # Plotting + im0 = axs[0].pcolormesh( + zax[::nr]/micron,xax[::nr]/micron,Ex.T[::nr,::nr], + vmin=-Exmax, vmax=Exmax, + cmap="RdBu", rasterized=True) + + plt.colorbar(im0,ax=axs[00], label=r"$E_x$ (V/m)") + + im1 = axs[1].pcolormesh( + zax[::nr]/micron,xax[::nr]/micron,By.T[::nr,::nr], + vmin=-Bymax, vmax=Bymax, + cmap="RdBu", rasterized=True) + plt.colorbar(im1,ax=axs[1], label=r"$B_y$ (T)") + + im2 = axs[2].pcolormesh( + zax[::nr]/micron,xax[::nr]/micron,Ez.T[::nr,::nr], + vmin=-Ezmax, vmax=Ezmax, + cmap="RdBu", rasterized=True) + plt.colorbar(im2,ax=axs[2],label=r"$E_z$ (V/m)") + + # Axis labels and title + for ax in axs: + ax.set_aspect(1.0) + ax.set_ylabel(r"$x$ ($\mu$m)") + for ax in axs[:-1]: + ax.set_xticklabels([]) + axs[2].set_xlabel(r"$z$ ($\mu$m)") + fig.suptitle(f"Iteration: {it}, Time: {time/1e-15:.1f} fs") + + plt.tight_layout() + + plt.savefig(f"{out_dir}/fields_{it:06d}.png") + +def visualize_particle_histogram_iteration(diag_name="histuH", species="hydrogen", iteration=1000, out_dir="./analysis"): + + it = iteration + + if species == "hydrogen": + # proton rest energy in eV + mc2 = sc.m_p/sc.electron_volt * sc.c**2 + elif species == "electron": + mc2 = sc.m_e/sc.electron_volt * sc.c**2 + else: + raise NotImplementedError("The only implemented presets for this analysis script are `electron` or `hydrogen`.") + + fs = 1.e-15 + MeV = 1.e6 + + df = pd.read_csv(f"./diags/reducedfiles/{diag_name}.txt",delimiter=r'\s+') + # the columns look like this: + # #[0]step() [1]time(s) [2]bin1=0.000220() [3]bin2=0.000660() [4]bin3=0.001100() + + # matches words, strings surrounded by " ' ", dots, minus signs and e for scientific notation in numbers + nested_list = [re.findall(r"[\w'\.]+",col) for col in df.columns] + + index = pd.MultiIndex.from_tuples(nested_list, names=('column#', 'name', 'bin value')) + + df.columns = (index) + + steps = df.values[:, 0].astype(int) + ii = np.where(steps == it)[0][0] + time = df.values[:, 1] + data = df.values[:, 2:] + edge_vals = np.array([float(row[2]) for row in df.columns[2:]]) + edges_MeV = (np.sqrt(edge_vals**2 + 1)-1) * mc2 / MeV + + time_fs = time / fs + + fig,ax = plt.subplots(1,1) + + ax.plot(edges_MeV, data[ii, :]) + ax.set_yscale("log") + ax.set_ylabel(r"d$N$/d$\mathcal{E}$ (arb. u.)") + ax.set_xlabel(r"$\mathcal{E}$ (MeV)") + + fig.suptitle(f"{species} - Iteration: {it}, Time: {time_fs[ii]:.1f} fs") + + plt.tight_layout() + plt.savefig(f"./{out_dir}/{diag_name}_{it:06d}.png") + + +if __name__ == "__main__": + + # Argument parsing + parser = argparse.ArgumentParser(description='Visualize Laser-Ion Accelerator Densities and Fields') + parser.add_argument('-d', '--diag_dir', type=str, default='./diags/diag1', help='Directory containing density and field diagnostics') + parser.add_argument('-i', '--iteration', type=int, default=None, help='Specific iteration to visualize') + parser.add_argument('-hn', '--histogram_name', type=str, default='histuH', help='Name of histogram diagnostic to visualize') + parser.add_argument('-hs', '--histogram_species', type=str, default='hydrogen', help='Particle species in the visualized histogram diagnostic') + args = parser.parse_args() + + # Create analysis directory + analysis_dir = 'analysis' + create_analysis_dir(analysis_dir) + + # Loading the time series + ts = OpenPMDTimeSeries(args.diag_dir) + + if args.iteration is not None: + visualize_density_iteration(ts, args.iteration, analysis_dir) + visualize_field_iteration(ts, args.iteration, analysis_dir) + visualize_particle_histogram_iteration(args.histogram_name, args.histogram_species, args.iteration, analysis_dir) + else: + for it in ts.iterations: + visualize_density_iteration(ts, it, analysis_dir) + visualize_field_iteration(ts, it, analysis_dir) + visualize_particle_histogram_iteration(args.histogram_name, args.histogram_species, it, analysis_dir) diff --git a/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json b/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json index 1678e794683..81850192121 100644 --- a/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json +++ b/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json @@ -1,33 +1,36 @@ { - "lev=0": { - "Bx": 0.0, - "By": 11411047.898351602, - "Bz": 0.0, - "Ex": 2031641076084948.5, - "Ey": 0.0, - "Ez": 334777106874158.8, - "jx": 1.6602050255725978e+19, - "jy": 0.0, - "jz": 9.545451466990307e+18, - "rho": 67237998613.86053, - "rho_electrons": 17449705826002.385, - "rho_hydrogen": 17441792654743.146 - }, - "hydrogen": { - "particle_momentum_x": 2.2575255298569512e-18, - "particle_momentum_z": 1.0773391029868316e-18, - "particle_orig_x": 0.007763427734375002, - "particle_orig_z": 0.0351975439453125, - "particle_position_x": 0.007763034484095496, - "particle_position_y": 0.035195761508116416, - "particle_weight": 2.6792730619156992e+17 - }, "electrons": { "particle_momentum_x": 3.8463518636930147e-19, "particle_momentum_y": 0.0, "particle_momentum_z": 1.6287398567676136e-18, "particle_position_x": 0.008139313907538123, - "particle_position_y": 0.0308605164193468, + "particle_position_y": 0.0, + "particle_position_z": 0.0308605164193468, "particle_weight": 2.647983436428149e+17 + }, + "hydrogen": { + "particle_momentum_x": 2.2575255298569516e-18, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.0773391029868316e-18, + "particle_origX": 0.007763427734375002, + "particle_origZ": 0.0351975439453125, + "particle_position_x": 0.007763034484095496, + "particle_position_y": 0.0, + "particle_position_z": 0.035195761508116416, + "particle_weight": 2.6792730619156992e+17 + }, + "lev=0": { + "Bx": 0.0, + "By": 11411047.898351604, + "Bz": 0.0, + "Ex": 2031641076084948.8, + "Ey": 0.0, + "Ez": 334777106874158.75, + "jx": 1.660205025572598e+19, + "jy": 0.0, + "jz": 9.545451466990307e+18, + "rho": 67237998613.860535, + "rho_electrons": 17449705826002.387, + "rho_hydrogen": 17441792654743.145 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/Python_LaserIonAcc2d.json b/Regression/Checksum/benchmarks_json/Python_LaserIonAcc2d.json new file mode 100644 index 00000000000..baaf29bec59 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/Python_LaserIonAcc2d.json @@ -0,0 +1,34 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 11406516.737324182, + "Bz": 0.0, + "Ex": 2032785026975781.0, + "Ey": 0.0, + "Ez": 315269171271122.75, + "jx": 1.6224680921022521e+19, + "jy": 0.0, + "jz": 8.875043590548824e+18, + "rho": 61329928310.45101, + "rho_electrons": 17450329551552.684, + "rho_hydrogen": 17441799909767.852 + }, + "electrons": { + "particle_position_x": 0.008152777855374778, + "particle_position_y": 0.0, + "particle_position_z": 0.030714485915215906, + "particle_momentum_x": 3.594599848069649e-19, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.6339243089017708e-18, + "particle_weight": 2.6507336926909222e+17 + }, + "hydrogen": { + "particle_position_x": 0.008197892199782453, + "particle_position_y": 0.0, + "particle_position_z": 0.0365646600930625, + "particle_momentum_x": 2.2464037026384752e-18, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.0873094324185116e-18, + "particle_weight": 2.703612070965676e+17 + } +} \ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index c9feea370af..4d1435e61de 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2184,7 +2184,8 @@ doVis = 0 [LaserIonAcc2d] buildDir = . inputFile = Examples/Physics_applications/laser_ion/inputs_2d -runtime_params = amr.n_cell=384 512 max_step=100 +outputFile = LaserIonAcc2d_plt +runtime_params = dim = 2 addToCompileString = USE_OPENPMD=TRUE cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_OPENPMD=ON @@ -2195,7 +2196,7 @@ useOMP = 1 numthreads = 1 compileTest = 0 doVis = 0 -analysisRoutine = Examples/analysis_default_regression.py +analysisRoutine = Examples/analysis_default_openpmd_regression.py [LaserOnFine] buildDir = . @@ -3393,6 +3394,25 @@ compareParticles = 1 particleTypes = electrons beam analysisRoutine = Examples/analysis_default_regression.py +[Python_LaserIonAcc2d] +buildDir = . +inputFile = Examples/Physics_applications/laser_ion/PICMI_inputs_2d.py +outputFile= diags/Python_LaserIonAcc2d_plt +runtime_params = +customRunCmd = python3 PICMI_inputs_2d.py +dim = 2 +addToCompileString = USE_OPENPMD=TRUE +cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_OPENPMD=ON -DWarpX_APP=OFF -DWarpX_PYTHON=ON +target = pip_install +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/analysis_default_openpmd_regression.py + [Python_LoadExternalField3D] buildDir = . inputFile = Examples/Tests/LoadExternalField/PICMI_inputs_3d.py