diff options
15 files changed, 724 insertions, 15 deletions
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 3bbe03722..8bd8c6cdf 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -1961,6 +1961,27 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a Default is ``<diag_name>.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz``. Note that the fields are averaged on the cell centers before they are written to file. +* ``<diag_name>.particle_fields_to_plot`` (list of `strings`, optional) + Names of per-cell averages of particle properties to calculate and output as additional fields. + Note that these averages do not respect the particle shape factor, and instead use nearest-grid point interpolation. + Default is none. + Parser functions for these field names are specified by ``<diag_name>.particle_fields.<field_name>(x,y,z,ux,uy,uz)``. + +* ``<diag_name>.particle_fields_species`` (list of `strings`, optional) + Species for which to calculate ``particle_fields_to_plot``. + Fields will be calculated separately for each specified species. + The default is a list of all of the available particle species. + +* ``<diag_name>.particle_fields.<field_name>(x,y,z,ux,uy,uz)`` (parser `string`) + Parser function to be calculated for each particle and averaged per cell. The field written is + + .. math:: + + \texttt{<field_name>_<species>} = \frac{\sum_{i=1}^N w_i \, f(x_i,y_i,z_i,u_{x,i},u_{y,i},u_{z,i})}{\sum_{i=1}^N w_i} + + where the sums are over all particles of type ``<species>`` in a cell (ignoring the particle shape factor), :math:`w_i` is the particle weight, :math:`f()` is the parser function, and :math:`(x_i,y_i,z_i)` are particle positions in units of a meter. + In 1D or 2D, the particle coordinates will follow the WarpX convention. :math:`(u_{x,i},u_{y,i},u_{z,i})` are components of the particle four-velocity. :math:`u = \gamma v/c`, :math:`\gamma` is the Lorentz factor, :math:`v` is the particle velocity, and :math:`c` is the speed of light. + * ``<diag_name>.plot_raw_fields`` (`0` or `1`) optional (default `0`) By default, the fields written in the plot files are averaged on the cell centers. When ``<diag_name>.plot_raw_fields = 1``, then the raw (i.e. non-averaged) diff --git a/Examples/Tests/particle_fields_diags/analysis_particle_diags.py b/Examples/Tests/particle_fields_diags/analysis_particle_diags.py new file mode 100755 index 000000000..3a77cbfc5 --- /dev/null +++ b/Examples/Tests/particle_fields_diags/analysis_particle_diags.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python3 + +# Copyright 2019-2022 Luca Fedeli, Yinjian Zhao, Hannah Klion +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +# This script tests the reduced particle diagnostics. +# The setup is a uniform plasma with electrons, protons and photons. +# Various particle and field quantities are written to file using the reduced diagnostics +# and compared with the corresponding quantities computed from the data in the plotfiles. + +import analysis_particle_diags_impl as an + +an.do_analysis(single_precision = False) diff --git a/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py b/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py new file mode 100755 index 000000000..e6dfc6a85 --- /dev/null +++ b/Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py @@ -0,0 +1,156 @@ +#!/usr/bin/env python3 + +# Copyright 2019-2022 Luca Fedeli, Yinjian Zhao, Hannah Klion +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +# This script tests the reduced particle diagnostics. +# The setup is a uniform plasma with electrons, protons and photons. +# Various particle and field quantities are written to file using the reduced diagnostics +# and compared with the corresponding quantities computed from the data in the plotfiles. + +import os +import sys + +import numpy as np +from scipy.constants import c +from scipy.constants import epsilon_0 as eps0 +from scipy.constants import m_e, m_p +from scipy.constants import mu_0 as mu0 +import yt + +sys.path.insert(1, '../../../../warpx/Regression/Checksum/') +import checksumAPI + + +def do_analysis(single_precision = False): + fn = sys.argv[1] + + ds = yt.load(fn) + ad = ds.all_data() + ad0 = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) + + #-------------------------------------------------------------------------------------------------- + # Part 1: get results from plotfiles (label '_yt') + #-------------------------------------------------------------------------------------------------- + + # Quantities computed from plotfiles + values_yt = dict() + + domain_size = ds.domain_right_edge.value - ds.domain_left_edge.value + dx = domain_size / ds.domain_dimensions + + # Electrons + x = ad['electrons', 'particle_position_x'].to_ndarray() + y = ad['electrons', 'particle_position_y'].to_ndarray() + z = ad['electrons', 'particle_position_z'].to_ndarray() + uz = ad['electrons', 'particle_momentum_z'].to_ndarray() / m_e / c + w = ad['electrons', 'particle_weight'].to_ndarray() + + x_ind = ((x - ds.domain_left_edge[0].value) / dx[0]).astype(int) + y_ind = ((y - ds.domain_left_edge[1].value) / dx[1]).astype(int) + z_ind = ((z - ds.domain_left_edge[2].value) / dx[2]).astype(int) + + zavg = np.zeros(ds.domain_dimensions) + uzavg = np.zeros(ds.domain_dimensions) + zuzavg = np.zeros(ds.domain_dimensions) + wavg = np.zeros(ds.domain_dimensions) + + for i_p in range(len(x)): + zavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += z[i_p] * w[i_p] + uzavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += uz[i_p] * w[i_p] + zuzavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += z[i_p] * uz[i_p] * w[i_p] + wavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += w[i_p] + + wavg_adj = np.where(wavg == 0, 1, wavg) + values_yt['electrons: zavg'] = zavg / wavg_adj + values_yt['electrons: uzavg'] = uzavg / wavg_adj + values_yt['electrons: zuzavg'] = zuzavg / wavg_adj + + # protons + x = ad['protons', 'particle_position_x'].to_ndarray() + y = ad['protons', 'particle_position_y'].to_ndarray() + z = ad['protons', 'particle_position_z'].to_ndarray() + uz = ad['protons', 'particle_momentum_z'].to_ndarray() / m_p / c + w = ad['protons', 'particle_weight'].to_ndarray() + + x_ind = ((x - ds.domain_left_edge[0].value) / dx[0]).astype(int) + y_ind = ((y - ds.domain_left_edge[1].value) / dx[1]).astype(int) + z_ind = ((z - ds.domain_left_edge[2].value) / dx[2]).astype(int) + + zavg = np.zeros(ds.domain_dimensions) + uzavg = np.zeros(ds.domain_dimensions) + zuzavg = np.zeros(ds.domain_dimensions) + wavg = np.zeros(ds.domain_dimensions) + + for i_p in range(len(x)): + zavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += z[i_p] * w[i_p] + uzavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += uz[i_p] * w[i_p] + zuzavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += z[i_p] * uz[i_p] * w[i_p] + wavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += w[i_p] + + wavg_adj = np.where(wavg == 0, 1, wavg) + values_yt['protons: zavg'] = zavg / wavg_adj + values_yt['protons: uzavg'] = uzavg / wavg_adj + values_yt['protons: zuzavg'] = zuzavg / wavg_adj + + # Photons (momentum in units of m_e c) + x = ad['photons', 'particle_position_x'].to_ndarray() + y = ad['photons', 'particle_position_y'].to_ndarray() + z = ad['photons', 'particle_position_z'].to_ndarray() + uz = ad['photons', 'particle_momentum_z'].to_ndarray() / m_e / c + w = ad['photons', 'particle_weight'].to_ndarray() + + x_ind = ((x - ds.domain_left_edge[0].value) / dx[0]).astype(int) + y_ind = ((y - ds.domain_left_edge[1].value) / dx[1]).astype(int) + z_ind = ((z - ds.domain_left_edge[2].value) / dx[2]).astype(int) + + zavg = np.zeros(ds.domain_dimensions) + uzavg = np.zeros(ds.domain_dimensions) + zuzavg = np.zeros(ds.domain_dimensions) + wavg = np.zeros(ds.domain_dimensions) + + for i_p in range(len(x)): + zavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += z[i_p] * w[i_p] + uzavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += uz[i_p] * w[i_p] + zuzavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += z[i_p] * uz[i_p] * w[i_p] + wavg[x_ind[i_p],y_ind[i_p],z_ind[i_p]] += w[i_p] + + wavg_adj = np.where(wavg == 0, 1, wavg) + values_yt['photons: zavg'] = zavg / wavg_adj + values_yt['photons: uzavg'] = uzavg / wavg_adj + values_yt['photons: zuzavg'] = zuzavg / wavg_adj + + + values_rd = dict() + # Load reduced particle diagnostic data from plotfiles + values_rd['electrons: zavg'] = ad0[('boxlib','z_electrons')] + values_rd['protons: zavg'] = ad0[('boxlib','z_protons')] + values_rd['photons: zavg'] = ad0[('boxlib','z_photons')] + + values_rd['electrons: uzavg'] = ad0[('boxlib','uz_electrons')] + values_rd['protons: uzavg'] = ad0[('boxlib','uz_protons')] + values_rd['photons: uzavg'] = ad0[('boxlib','uz_photons')] + + values_rd['electrons: zuzavg'] = ad0[('boxlib','zuz_electrons')] + values_rd['protons: zuzavg'] = ad0[('boxlib','zuz_protons')] + values_rd['photons: zuzavg'] = ad0[('boxlib','zuz_photons')] + + #-------------------------------------------------------------------------------------------------- + # Part 3: compare values from plotfiles and diagnostics and print output + #-------------------------------------------------------------------------------------------------- + + error = dict() + tolerance = 5e-3 if single_precision else 1e-12 + + for k in values_yt.keys(): + # check that the zeros line up, since we'll be ignoring them in the error calculation + assert(np.all((values_yt[k] == 0) == (values_rd[k] == 0))) + error[k] = np.max(abs(values_yt[k] - values_rd[k])[values_yt[k] != 0] / abs(values_yt[k])[values_yt[k] != 0]) + print(k, 'relative error = ', error[k]) + assert(error[k] < tolerance) + + test_name = os.path.split(os.getcwd())[1] + checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Tests/particle_fields_diags/analysis_particle_diags_single.py b/Examples/Tests/particle_fields_diags/analysis_particle_diags_single.py new file mode 100755 index 000000000..56d98831e --- /dev/null +++ b/Examples/Tests/particle_fields_diags/analysis_particle_diags_single.py @@ -0,0 +1,16 @@ +#!/usr/bin/env python3 + +# Copyright 2019-2022 Luca Fedeli, Yinjian Zhao, Hannah Klion +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL + +# This script tests the reduced particle diagnostics. +# The setup is a uniform plasma with electrons, protons and photons. +# Various particle and field quantities are written to file using the reduced diagnostics +# and compared with the corresponding quantities computed from the data in the plotfiles. + +import analysis_particle_diags_impl as an + +an.do_analysis(single_precision = True) diff --git a/Examples/Tests/particle_fields_diags/inputs b/Examples/Tests/particle_fields_diags/inputs new file mode 100644 index 000000000..5684908df --- /dev/null +++ b/Examples/Tests/particle_fields_diags/inputs @@ -0,0 +1,78 @@ +# Maximum number of time steps +max_step = 200 + +# 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 = 32 + +# Maximum level in hierarchy +amr.max_level = 0 + +# Geometry +geometry.dims = 3 +geometry.prob_lo = -1. -1. -1. # physical domain +geometry.prob_hi = 1. 1. 1. + +# Boundary condition +boundary.field_lo = periodic periodic periodic +boundary.field_hi = periodic periodic periodic + +# Algorithms +algo.current_deposition = esirkepov +algo.field_gathering = energy-conserving # or momentum-conserving +warpx.use_filter = 1 +algo.maxwell_solver = yee # or ckc + +# Order of particle shape factors +algo.particle_shape = 1 + +# CFL +warpx.cfl = 0.99999 + +# Particles +particles.species_names = electrons protons photons +particles.photon_species = photons + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = 1 1 1 +electrons.profile = constant +electrons.density = 1.e14 # number of electrons per m^3 +electrons.momentum_distribution_type = gaussian +electrons.ux_th = 0.035 +electrons.uy_th = 0.035 +electrons.uz_th = 0.035 + +protons.charge = q_e +protons.mass = m_p +protons.injection_style = "NUniformPerCell" +protons.num_particles_per_cell_each_dim = 1 1 1 +protons.profile = constant +protons.density = 1.e14 # number of protons per m^3 +protons.momentum_distribution_type = at_rest + +photons.species_type = "photon" +photons.injection_style = "NUniformPerCell" +photons.num_particles_per_cell_each_dim = 1 1 1 +photons.profile = constant +photons.density = 1.e14 # number of protons per m^3 +photons.momentum_distribution_type = gaussian +photons.ux_th = 0.2 +photons.uy_th = 0.2 +photons.uz_th = 0.2 + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 200 +diag1.diag_type = Full +diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_protons + +diag1.particle_fields_to_plot = z uz zuz +diag1.particle_fields_species = electrons protons photons +diag1.particle_fields.z(x,y,z,ux,uy,uz) = z +diag1.particle_fields.uz(x,y,z,ux,uy,uz) = uz +diag1.particle_fields.zuz(x,y,z,ux,uy,uz) = z * uz diff --git a/Regression/Checksum/benchmarks_json/particle_fields_diags.json b/Regression/Checksum/benchmarks_json/particle_fields_diags.json new file mode 100644 index 000000000..4629b0b27 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/particle_fields_diags.json @@ -0,0 +1,58 @@ +{ + "electrons": { + "particle_cpu": 16384.0, + "particle_id": 268451840.0, + "particle_momentum_x": 2.4335130953142086e-19, + "particle_momentum_y": 2.463314849025226e-19, + "particle_momentum_z": 2.4452967447264526e-19, + "particle_position_x": 16386.79272675649, + "particle_position_y": 16383.137717233834, + "particle_position_z": 16385.771013436024, + "particle_weight": 800000000000000.0 + }, + "lev=0": { + "Bx": 0.08405082842287068, + "By": 0.08395442339587296, + "Bz": 0.08318206628870235, + "Ex": 102195850.94145432, + "Ey": 106377257.8660376, + "Ez": 102627869.95880052, + "jx": 714393.4493262022, + "jy": 739611.1829573747, + "jz": 719566.2651192, + "rho": 0.02721945765330138, + "rho_electrons": 0.5250012394291199, + "rho_protons": 0.5250012394291199, + "uz_electrons": 502.0697557311614, + "uz_photons": 2868.9771458584096, + "uz_protons": 0.27498556151402315, + "z_electrons": 10620.91075174126, + "z_photons": 10315.870910754074, + "z_protons": 16383.994708420258, + "zuz_electrons": 251.42139943815852, + "zuz_photons": 1423.0530482617983, + "zuz_protons": 0.13812948570246372 + }, + "photons": { + "particle_cpu": 16384.0, + "particle_id": 1342193664.0, + "particle_momentum_x": 1.43264587339745e-18, + "particle_momentum_y": 1.420927243614407e-18, + "particle_momentum_z": 1.4314382225453842e-18, + "particle_position_x": 16274.031402786917, + "particle_position_y": 16374.776556959983, + "particle_position_z": 16308.352100156182, + "particle_weight": 800000000000000.0 + }, + "protons": { + "particle_cpu": 16384.0, + "particle_id": 805322752.0, + "particle_momentum_x": 1.4305311394194743e-19, + "particle_momentum_y": 1.4342178041433115e-19, + "particle_momentum_z": 1.378886053708302e-19, + "particle_position_x": 16384.020927263282, + "particle_position_y": 16384.01049132875, + "particle_position_z": 16383.994708420258, + "particle_weight": 800000000000000.0 + } +}
\ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/particle_fields_diags_single_precision.json b/Regression/Checksum/benchmarks_json/particle_fields_diags_single_precision.json new file mode 100644 index 000000000..bfbafee9e --- /dev/null +++ b/Regression/Checksum/benchmarks_json/particle_fields_diags_single_precision.json @@ -0,0 +1,58 @@ +{ + "electrons": { + "particle_cpu": 16384.0, + "particle_id": 268451840.0, + "particle_momentum_x": 2.441720970290429e-19, + "particle_momentum_y": 2.4564842258504588e-19, + "particle_momentum_z": 2.434157527179439e-19, + "particle_position_x": 16379.247611900064, + "particle_position_y": 16382.836088184762, + "particle_position_z": 16382.833910756759, + "particle_weight": 800000003014656.0 + }, + "lev=0": { + "Bx": 0.08424076243170475, + "By": 0.08475179390317322, + "Bz": 0.08796453154337058, + "Ex": 106872501.50568771, + "Ey": 107805476.44736862, + "Ez": 107783029.62442589, + "jx": 728220.0582114309, + "jy": 736403.7788231075, + "jz": 734813.6176411062, + "rho": 0.027750853105895423, + "rho_electrons": 0.5250012279730072, + "rho_protons": 0.5250012273099856, + "uz_electrons": 499.5564858153293, + "uz_photons": 2823.3768876476906, + "uz_protons": 0.27726521408384874, + "z_electrons": 10663.295547661037, + "z_photons": 10320.570131806278, + "z_protons": 16384.033098796383, + "zuz_electrons": 249.21871114507326, + "zuz_photons": 1421.2299280466445, + "zuz_protons": 0.1413706733222889 + }, + "photons": { + "particle_cpu": 16384.0, + "particle_id": 1342193664.0, + "particle_momentum_x": 1.428291590249666e-18, + "particle_momentum_y": 1.4222174024686332e-18, + "particle_momentum_z": 1.4246585206989104e-18, + "particle_position_x": 16376.414849636145, + "particle_position_y": 16409.713781019993, + "particle_position_z": 16378.407608640264, + "particle_weight": 800000003014656.0 + }, + "protons": { + "particle_cpu": 16384.0, + "particle_id": 805322752.0, + "particle_momentum_x": 1.4104800123456525e-19, + "particle_momentum_y": 1.4120344417599842e-19, + "particle_momentum_z": 1.3903171547264748e-19, + "particle_position_x": 16383.951009677723, + "particle_position_y": 16383.991230150685, + "particle_position_z": 16384.03309544921, + "particle_weight": 800000003014656.0 + } +}
\ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index d941aa225..f129516c0 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -2196,6 +2196,42 @@ doVis = 0 compareParticles = 0 analysisRoutine = Examples/Tests/reduced_diags/analysis_reduced_diags_loadbalancecosts.py +[particle_fields_diags] +buildDir = . +inputFile = Examples/Tests/particle_fields_diags/inputs +aux1File = Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 +dim = 3 +addToCompileString = +cmakeSetupOpts = -DWarpX_DIMS=3 +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Tests/particle_fields_diags/analysis_particle_diags.py + +[particle_fields_diags_single_precision] +buildDir = . +inputFile = Examples/Tests/particle_fields_diags/inputs +aux1File = Examples/Tests/particle_fields_diags/analysis_particle_diags_impl.py +runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_initial_conditions=1 +dim = 3 +addToCompileString = PRECISION=FLOAT USE_SINGLE_PRECISION_PARTICLES=TRUE +cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_PRECISION=SINGLE +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Tests/particle_fields_diags/analysis_particle_diags_single.py + [galilean_2d_psatd] buildDir = . inputFile = Examples/Tests/galilean/inputs_2d diff --git a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt index 5e085db4d..b2d0c096d 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt +++ b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt @@ -8,4 +8,5 @@ target_sources(WarpX PartPerGridFunctor.cpp BackTransformFunctor.cpp BackTransformParticleFunctor.cpp + ParticleReductionFunctor.cpp ) diff --git a/Source/Diagnostics/ComputeDiagFunctors/Make.package b/Source/Diagnostics/ComputeDiagFunctors/Make.package index 7a07273e2..611ccb60d 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/Make.package +++ b/Source/Diagnostics/ComputeDiagFunctors/Make.package @@ -6,5 +6,6 @@ CEXE_sources += DivEFunctor.cpp CEXE_sources += RhoFunctor.cpp CEXE_sources += BackTransformFunctor.cpp CEXE_sources += BackTransformParticleFunctor.cpp +CEXE_sources += ParticleReductionFunctor.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ComputeDiagFunctors diff --git a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H new file mode 100644 index 000000000..81172a423 --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H @@ -0,0 +1,50 @@ +#ifndef WARPX_PARTICLEREDUCTIONFUNCTOR_H_ +#define WARPX_PARTICLEREDUCTIONFUNCTOR_H_ + +#include "ComputeDiagFunctor.H" + +#include <AMReX_Parser.H> + +#include <AMReX_BaseFwd.H> + +#include <memory> +#include <string> + +/** + * \brief Functor to calculate per-cell averages of particle properties. + */ +class +ParticleReductionFunctor final : public ComputeDiagFunctor +{ +public: + /** Constructor. + * \param[in] mf_src source multifab. Must be nullptr as no source MF is needed + * to compute the number of particles per cell, banane. + * \param[in] lev level of multifab. + * \param[in] crse_ratio for interpolating field values from simulation MultiFabs + to the output diagnostic MultiFab mf_dst. + * \param[in] fn_str parser string that describes the function to apply to particles + * and then average over each cell + * \param[in] ispec index of the species over which to calculate the average + * \param[in] ncomp Number of component of mf_src to cell-center in dst multifab. + */ + ParticleReductionFunctor(const amrex::MultiFab * const mf_src, const int lev, + const amrex::IntVect crse_ratio, const std::string fn_str, const int ispec, const int ncomp=1); + + /** \brief Compute the average of the function m_map_fn over each grid cell. + * + * \param[out] mf_dst output MultiFab where the result is written + * \param[in] dcomp first component of mf_dst in which cell-centered + * data is stored + */ + virtual void operator()(amrex::MultiFab& mf_dst, const int dcomp, const int /*i_buffer=0*/) const override; +private: + int const m_lev; /**< level on which mf_src is defined */ + int const m_ispec; /**< index of species to average over */ + /** Parser function to be averaged by the functor. Arguments: x, y, z, ux, uy, uz */ + std::unique_ptr<amrex::Parser> m_map_fn_parser; + /** Compiled #m_map_fn_parser */ + amrex::ParserExecutor<6> m_map_fn; +}; + +#endif // WARPX_PARTICLEREDUCTIONFUNCTOR_H_ diff --git a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp new file mode 100644 index 000000000..c5f9028f5 --- /dev/null +++ b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp @@ -0,0 +1,130 @@ + +#include "ParticleReductionFunctor.H" + +#include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H" +#include "Particles/MultiParticleContainer.H" +#include "Particles/WarpXParticleContainer.H" +#include "Utils/CoarsenIO.H" +#include "WarpX.H" + +#include <AMReX_Array.H> +#include <AMReX_BLassert.H> +#include <AMReX_IntVect.H> +#include <AMReX_MultiFab.H> +#include <AMReX_REAL.H> + +using namespace amrex::literals; + +ParticleReductionFunctor::ParticleReductionFunctor (const amrex::MultiFab* mf_src, const int lev, + const amrex::IntVect crse_ratio, const std::string fn_str, const int ispec, const int ncomp) + : ComputeDiagFunctor(ncomp, crse_ratio), m_lev(lev), m_ispec(ispec) +{ + // mf_src will not be used, let's make sure it's null. + AMREX_ALWAYS_ASSERT(mf_src == nullptr); + // Write only in one output component. + AMREX_ALWAYS_ASSERT(ncomp == 1); + + // Allocate and compile a parser based on the input string fn_str + m_map_fn_parser = std::make_unique<amrex::Parser>(makeParser( + fn_str, {"x", "y", "z", "ux", "uy", "uz"})); + m_map_fn = m_map_fn_parser->compile<6>(); +} + +void +ParticleReductionFunctor::operator() (amrex::MultiFab& mf_dst, const int dcomp, const int /*i_buffer*/) const +{ + auto& warpx = WarpX::GetInstance(); + // Guard cell is set to 1 for generality. However, for a cell-centered + // output Multifab, mf_dst, the guard-cell data is not needed especially considering + // the operations performend in the CoarsenAndInterpolate function. + constexpr int ng = 1; + // Temporary cell-centered, multi-component MultiFab for storing particles per cell. + amrex::MultiFab red_mf(warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), 1, ng); + amrex::MultiFab ppc_mf(warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), 1, ng); + // Set value to 0, and increment the value in each cell with ppc. + red_mf.setVal(0._rt); + ppc_mf.setVal(0._rt); + auto& pc = warpx.GetPartContainer().GetParticleContainer(m_ispec); + // Copy over compiled parser function so that it can be captured by value in the lambda + auto map_fn = m_map_fn; + ParticleToMesh(pc, red_mf, m_lev, + [=] AMREX_GPU_DEVICE (const WarpXParticleContainer::SuperParticleType& p, + amrex::Array4<amrex::Real> const& out_array, + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo, + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi) + { + // Get position in WarpX convention to use in parser. Will be different from + // p.pos() for 1D and 2D simulations. + amrex::ParticleReal xw = 0._rt, yw = 0._rt, zw = 0._rt; + get_particle_position(p, xw, yw, zw); + + // Get position in AMReX convention to calculate corresponding index. + // Ideally this will be replaced with the AMReX NGP interpolator + // Always do x direction. No RZ case because it's not implemented, and code + // will have aborted + int ii = 0, jj = 0, kk = 0; + amrex::ParticleReal x = p.pos(0); + amrex::Real lx = (x - plo[0]) * dxi[0]; + ii = static_cast<int>(amrex::Math::floor(lx)); +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_3D) + amrex::ParticleReal y = p.pos(1); + amrex::Real ly = (y - plo[1]) * dxi[1]; + jj = static_cast<int>(amrex::Math::floor(ly)); +#endif +#if defined(WARPX_DIM_3D) + amrex::ParticleReal z = p.pos(2); + amrex::Real lz = (z - plo[2]) * dxi[2]; + kk = static_cast<int>(amrex::Math::floor(lz)); +#endif + + // Fix dimensions since parser assumes u = gamma * v / c + amrex::ParticleReal ux = p.rdata(PIdx::ux) / PhysConst::c; + amrex::ParticleReal uy = p.rdata(PIdx::uy) / PhysConst::c; + amrex::ParticleReal uz = p.rdata(PIdx::uz) / PhysConst::c; + amrex::Real value = map_fn(xw, yw, zw, ux, uy, uz); + amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), p.rdata(PIdx::w) * value); + }); + // Add the weight for each particle -- total number of particles of this species + ParticleToMesh(pc, ppc_mf, m_lev, + [=] AMREX_GPU_DEVICE (const WarpXParticleContainer::SuperParticleType& p, + amrex::Array4<amrex::Real> const& out_array, + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo, + amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi) + { + + // Get position in AMReX convention to calculate corresponding index. + // Ideally this will be replaced with the AMReX NGP interpolator + // Always do x direction. No RZ case because it's not implemented, and code + // will have aborted + int ii = 0, jj = 0, kk = 0; + amrex::ParticleReal x = p.pos(0); + amrex::Real lx = (x - plo[0]) * dxi[0]; + ii = static_cast<int>(amrex::Math::floor(lx)); +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_3D) + amrex::ParticleReal y = p.pos(1); + amrex::Real ly = (y - plo[1]) * dxi[1]; + jj = static_cast<int>(amrex::Math::floor(ly)); +#endif +#if defined(WARPX_DIM_3D) + amrex::ParticleReal z = p.pos(2); + amrex::Real lz = (z - plo[2]) * dxi[2]; + kk = static_cast<int>(amrex::Math::floor(lz)); +#endif + amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), p.rdata(PIdx::w)); + }); + // Divide value by number of particles for average. Set average to zero if there are no particles + for (amrex::MFIter mfi(red_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) + { + const amrex::Box& box = mfi.tilebox(); + amrex::Array4<amrex::Real> const& a_red = red_mf.array(mfi); + amrex::Array4<amrex::Real const> const& a_ppc = ppc_mf.const_array(mfi); + amrex::ParallelFor(box, + [=] AMREX_GPU_DEVICE (int i, int j, int k) { + if (a_ppc(i,j,k,0) == 0) a_red(i,j,k,0) = 0; + else a_red(i,j,k,0) = a_red(i,j,k,0) / a_ppc(i,j,k,0); + }); + } + + // Coarsen and interpolate from ppc_mf to the output diagnostic MultiFab, mf_dst. + CoarsenIO::Coarsen(mf_dst, red_mf, dcomp, 0, nComp(), 0, m_crse_ratio); +} diff --git a/Source/Diagnostics/Diagnostics.H b/Source/Diagnostics/Diagnostics.H index a88e6221c..4b90cdb22 100644 --- a/Source/Diagnostics/Diagnostics.H +++ b/Source/Diagnostics/Diagnostics.H @@ -134,9 +134,20 @@ protected: /** Index of diagnostics in MultiDiagnostics::alldiags */ int m_diag_index; /** Names of each component requested by the user. - * in cylindrical geometry, this list is appended with - * automatically-constructed names for all modes of all fields */ + * The list is appended with the average particle fields, if used. + * In cylindrical geometry, this list is also appended with + * automatically-constructed names for all modes of all fields.*/ amrex::Vector< std::string > m_varnames; + /** Names of plotfile fields requested by the user */ + amrex::Vector< std::string > m_varnames_fields; + /** Names of averaged particle properties to output */ + amrex::Vector< std::string > m_pfield_varnames; + /** Names of species for which to output average particle fields */ + std::vector< std::string > m_pfield_species; + /** Species indices corresponding to elements of m_pfield_varnames */ + std::vector< int > m_pfield_species_index; + /** Map connecting particle field variable names to the corresponding parser strings */ + std::map< std::string, std::string > m_pfield_strings = {}; /** If true, a dump is performed at the last timestep regardless of the required dump timesteps */ bool m_dump_last_timestep = true; /** format for output files, "plotfile" or "openpmd" or "sensei" or "ascent" diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index fb41bebf1..b36011c1c 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/Source/Diagnostics/Diagnostics.cpp @@ -56,27 +56,27 @@ Diagnostics::BaseReadParameters () pp_diag_name.query("dump_last_timestep", m_dump_last_timestep); // Query list of grid fields to write to output - bool varnames_specified = pp_diag_name.queryarr("fields_to_plot", m_varnames); + bool varnames_specified = pp_diag_name.queryarr("fields_to_plot", m_varnames_fields); if (!varnames_specified){ - m_varnames = {"Ex", "Ey", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz"}; + m_varnames_fields = {"Ex", "Ey", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz"}; } // Sanity check if user requests to plot phi - if (WarpXUtilStr::is_in(m_varnames, "phi")){ + if (WarpXUtilStr::is_in(m_varnames_fields, "phi")){ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( warpx.do_electrostatic==ElectrostaticSolverAlgo::LabFrame, "plot phi only works if do_electrostatic = labframe"); } // Sanity check if user requests to plot F - if (WarpXUtilStr::is_in(m_varnames, "F")){ + if (WarpXUtilStr::is_in(m_varnames_fields, "F")){ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( warpx.do_dive_cleaning, "plot F only works if warpx.do_dive_cleaning = 1"); } // G can be written to file only if WarpX::do_divb_cleaning = 1 - if (WarpXUtilStr::is_in(m_varnames, "G")) + if (WarpXUtilStr::is_in(m_varnames_fields, "G")) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( warpx.do_divb_cleaning, "G can be written to file only if warpx.do_divb_cleaning = 1"); @@ -85,9 +85,74 @@ Diagnostics::BaseReadParameters () // If user requests to plot proc_number for a serial run, // delete proc_number from fields_to_plot if (amrex::ParallelDescriptor::NProcs() == 1){ - m_varnames.erase( - std::remove(m_varnames.begin(), m_varnames.end(), "proc_number"), - m_varnames.end()); + m_varnames_fields.erase( + std::remove(m_varnames_fields.begin(), m_varnames_fields.end(), "proc_number"), + m_varnames_fields.end()); + } + + // Get names of particle quantities to average at each grid point + const bool pfield_varnames_specified = pp_diag_name.queryarr("particle_fields_to_plot", m_pfield_varnames); + if (!pfield_varnames_specified){ + m_pfield_varnames = {}; + } +#ifdef WARPX_DIM_RZ + if (m_pfield_varnames.size() != 0) { + amrex::Abort("Input error: cannot use particle_fields_to_plot not implemented for RZ"); + } +#endif + + + // Get parser strings for particle fields and generate map of parsers + std::string parser_str; + amrex::ParmParse pp_diag_pfield(m_diag_name + ".particle_fields"); + for (const auto& var : m_pfield_varnames) { + Store_parserString(pp_diag_pfield, (var + "(x,y,z,ux,uy,uz)").c_str(), parser_str); + if (parser_str != "") { + m_pfield_strings.insert({var, parser_str}); + } + else { + amrex::Abort("Input error: cannot find parser string for " + var + "." + + m_diag_name + ".particle_fields." + var + " in file"); + } + } + + // Names of all species in the simulation + m_all_species_names = warpx.GetPartContainer().GetSpeciesNames(); + + // Get names of species to average at each grid point + const bool pfield_species_specified = pp_diag_name.queryarr("particle_fields_species", m_pfield_species); + if (!pfield_species_specified){ + m_pfield_species = m_all_species_names; + } + + // Check that species names specified in m_pfield_species are valid + bool p_species_name_is_wrong; + // Loop over all species specified above + for (const auto& species : m_pfield_species) { + // Boolean used to check if species name was misspelled + p_species_name_is_wrong = true; + // Loop over all species + for (int i = 0, n = int(m_all_species_names.size()); i < n; i++) { + if (species == m_all_species_names[i]) { + // Store species index: will be used in ParticleReductionFunctor to calculate + // averages for this species + m_pfield_species_index.push_back(i); + p_species_name_is_wrong = false; + } + } + // If species name was misspelled, abort with error message + if (p_species_name_is_wrong) { + amrex::Abort("Input error: string " + species + " in " + m_diag_name + + ".particle_fields_species does not match any species"); + } + } + + m_varnames = m_varnames_fields; + // Generate names of averaged particle fields and append to m_varnames + for (int ivar=0; ivar<m_pfield_varnames.size(); ivar++) { + for (int ispec=0; ispec < int(m_pfield_species.size()); ispec++) { + m_varnames.push_back(m_pfield_varnames[ivar] + '_' + m_pfield_species[ispec]); + } } // Read user-defined physical extents for the output and store in m_lo and m_hi. @@ -139,8 +204,6 @@ Diagnostics::BaseReadParameters () // Names of species to write to output bool species_specified = pp_diag_name.queryarr("species", m_output_species_names); - // Names of all species in the simulation - m_all_species_names = warpx.GetPartContainer().GetSpeciesNames(); // Auxiliary variables std::string species; @@ -175,6 +238,8 @@ Diagnostics::BaseReadParameters () bool checkpoint_compatibility = false; if (m_format == "checkpoint"){ if ( varnames_specified == false && + pfield_varnames_specified == false && + pfield_species_specified == false && lo_specified == false && hi_specified == false && cr_specified == false && diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index fb7a055f6..79791ae8e 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/Source/Diagnostics/FullDiagnostics.cpp @@ -5,6 +5,7 @@ #include "ComputeDiagFunctors/DivEFunctor.H" #include "ComputeDiagFunctors/PartPerCellFunctor.H" #include "ComputeDiagFunctors/PartPerGridFunctor.H" +#include "ComputeDiagFunctors/ParticleReductionFunctor.H" #include "ComputeDiagFunctors/RhoFunctor.H" #include "Diagnostics/Diagnostics.H" #include "Diagnostics/ParticleDiag/ParticleDiag.H" @@ -408,10 +409,13 @@ FullDiagnostics::InitializeFieldFunctors (int lev) // Species index to loop over species that dump rho per species int i = 0; - m_all_field_functors[lev].resize( m_varnames.size() ); + const auto nvar = static_cast<int>(m_varnames_fields.size()); + const auto nspec = static_cast<int>(m_pfield_species.size()); + const auto ntot = static_cast<int>(nvar + m_pfield_varnames.size() * nspec); + + m_all_field_functors[lev].resize(ntot); // Fill vector of functors for all components except individual cylindrical modes. - const auto n = static_cast<int>(m_all_field_functors[lev].size()); - for (int comp=0; comp<n; comp++){ + for (int comp=0; comp<nvar; comp++){ if ( m_varnames[comp] == "Ex" ){ m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 0), lev, m_crse_ratio); } else if ( m_varnames[comp] == "Ey" ){ @@ -456,6 +460,14 @@ FullDiagnostics::InitializeFieldFunctors (int lev) amrex::Abort("Error: " + m_varnames[comp] + " is not a known field output type"); } } + // Add functors for average particle data for each species + for (int pcomp=0; pcomp<int(m_pfield_varnames.size()); pcomp++) { + std::string varname = m_pfield_varnames[pcomp]; + for (int ispec=0; ispec<int(m_pfield_species.size()); ispec++) { + m_all_field_functors[lev][nvar + pcomp * nspec + ispec] = std::make_unique<ParticleReductionFunctor>(nullptr, + lev, m_crse_ratio, m_pfield_strings[varname], m_pfield_species_index[ispec]); + } + } AddRZModesToDiags( lev ); } |