Add Plasma Histogram Diagnostic

CMakeLists.txt

@@ -254,6 +254,7 @@ if(BUILD_TESTING)
             linear_wake.SI.1Rank
             gaussian_linear_wake.normalized.1Rank
             gaussian_linear_wake.SI.1Rank
+            histogram.2Rank
             reset.2Rank
             beam_in_vacuum.SI.1Rank
             beam_in_vacuum.normalized.1Rank

docs/source/run/parameters.rst

@@ -1113,7 +1113,7 @@ Field diagnostics
 
 * ``<diag name> or diagnostic.base_geometry`` (`string`) optional (default `level_0`)
     Which geometry the diagnostics should be based on.
-    Available geometries are `level_0`, `level_1`, `level_2` and `laser`,
+    Available geometries are `level_0`, `level_1`, `level_2`, `laser` and `histogram`,
     depending on if MR or a laser is used.
     If ``<diag name>`` is equal to ``lev0 lev1 lev2 laser_diag``, the default for this parameter
     becomes ``level_0 level_1 level_2 laser`` respectively.
@@ -1192,6 +1192,56 @@ Field diagnostics
     this option specifies the shape order of the deposited fields.
     Currently, 0,1,2,3 are implemented.
 
+Particle Histogram diagnostics
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+This diagnostic allows the direct computation of histograms of arbitrary particle
+quantities during the simulation runtime. It supports both plasma and beam particles.
+It is part of the standard field diagnostic and is enabled by setting
+``<diag name>.base_geometry = histogram``. Histograms may have one or two user-defined axes.
+Optionally, the simulation z-axis can be included as an additional axis or integrated over.
+All field diagnostic parameters apply, except ``field_data`` and ``diag_type``.
+For example, ``patch_lo`` and ``patch_hi`` can be used to restrict the particles included in
+the histogram in coordinate space, but do not modify the histogram axes.
+Each particle species produces a separate histogram that is output as
+``<species name>_<diag name>``. No unit conversion or normalization by cell volume is applied.
+Particles outside the histogram bounds are discarded.
+
+* ``<diag name>.hist_species_names`` (`string`)
+    List of species to include. Can be beam and/or plasma species.
+
+* ``<diag name>.hist_num_bins`` (`int` or 2 `int`)
+    Number of bins per histogram axis.
+
+* ``<diag name>.hist_bins_lo`` (`flaot` or 2 `flaot`)
+    Lower bound of each histogram axis.
+
+* ``<diag name>.hist_bins_hi`` (`flaot` or 2 `flaot`)
+    Upper bound of each histogram axis.
+
+* ``<diag name>.hist_function`` (`string`)
+    Parser expression defining the first histogram axis as a function of particle properties:
+    ``x``, ``y``, ``z``, ``ux``, ``uy``, ``uz``, ``ga_psi``, ``w``, ``ion_lev``.
+    Here, ``ga_psi`` is the quasi-static weighting factor for plasma particles,
+    ``ion_lev`` is the ionization level (for ionizable species), and ``w`` is the
+    macro-particle weight.
+
+* ``<diag name>.hist_function2`` (`string`) optional
+    Parser expression defining the second histogram axis. Uses the same variables as
+    ``<diag name>.hist_function``.
+
+* ``<diag name>.hist_weight`` (`string`) optional (default `w`)
+    Parser expression defining the weight contributed by each particle to the histogram.
+    For plasma particles, this weight is automatically multiplied by ``ga_psi``
+    to obtain the physical particle weight. Uses the same variables as
+    ``<diag name>.hist_function``. This can also be used to filter particles.
+
+* ``<diag name>.hist_add_z_axis`` (`bool`) optional (default `false`)
+    Add the zeta axis from the simulation to the histogram output.
+    This is more efficient than adding z as a custom histogram axis using
+    ``hist_function`` or ``hist_function2``. If disabled the histogram contains data from
+    all z slices in the range given by ``patch_lo`` and ``patch_hi``.
+
 In-situ diagnostics
 ^^^^^^^^^^^^^^^^^^^
 

examples/blowout_wake/analysis_histogram.py

@@ -0,0 +1,24 @@
+#! /usr/bin/env python3
+
+# Copyright 2026
+#
+# This file is part of HiPACE++.
+#
+# Authors: AlexanderSinn
+# License: BSD-3-Clause-LBNL
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+
+
+ts = OpenPMDTimeSeries('histogram.2Rank')
+
+hist_uz1 = ts.get_field(field="beam_hist_z_uz1", iteration=ts.iterations[-1])[0]
+hist_uz2 = ts.get_field(field="beam_hist_z_uz2", iteration=ts.iterations[-1])[0]
+
+error = np.max(np.abs(hist_uz1 - hist_uz2.T)) / np.max(np.abs(hist_uz1))
+print("error =", error)
+assert(error < 1.e-20)
+
+del ts
+

examples/blowout_wake/inputs_histogram_normalized

@@ -0,0 +1,72 @@
+amr.n_cell = 64 64 100
+
+hipace.normalized_units=1
+
+amr.max_level = 0
+
+max_step = 3
+hipace.dt = 1
+
+hipace.depos_order_xy = 2
+
+boundary.field = Dirichlet
+boundary.particle = Periodic
+geometry.prob_lo = -8 -8 -6 # physical domain
+geometry.prob_hi =  8  8  6
+
+beams.names = beam
+beam.injection_type = fixed_weight_pdf
+beam.pdf = "exp(-0.5*(z/1.41)^2)"
+beam.density = 3.
+beam.u_mean = 0. 0. 2000
+beam.u_std = 0. 0. 0.
+beam.position_mean = 0 0
+beam.position_std = 0.3 0.3
+beam.num_particles = 1e5
+
+plasmas.names = plasma
+plasma.density(x,y,z) = 1.
+plasma.ppc = 2 2
+plasma.element = electron
+plasma.temperature_in_ev = 1
+
+diagnostic.diag_type = xz
+diagnostic.output_period = "current_step == max_step"
+diagnostic.beam_output_period = 0
+
+diagnostic.names = lev0 hist_z_uz1 hist_z_uz2 hist_temperature hist_all
+
+lev0.field_data = By ExmBy Ez Psi Sx chi jx_beam jz_beam rhomjz
+
+hist_z_uz1.base_geometry = histogram
+hist_z_uz1.hist_species_names = beam
+hist_z_uz1.hist_num_bins = 100 200
+hist_z_uz1.hist_bins_lo = -6 1980
+hist_z_uz1.hist_bins_hi =  6 2020
+hist_z_uz1.hist_function = "z"
+hist_z_uz1.hist_function2 = "uz"
+
+hist_z_uz2.base_geometry = histogram
+hist_z_uz2.hist_species_names = beam
+hist_z_uz2.hist_num_bins = 200
+hist_z_uz2.hist_bins_lo = 1980
+hist_z_uz2.hist_bins_hi = 2020
+hist_z_uz2.hist_function = "uz"
+hist_z_uz2.hist_add_z_axis = true
+
+hist_temperature.base_geometry = histogram
+hist_temperature.hist_species_names = plasma
+hist_temperature.hist_num_bins = 200
+hist_temperature.hist_bins_lo = 0
+hist_temperature.hist_bins_hi = 0.02
+hist_temperature.hist_function = "sqrt(ux^2 + uy^2 + uz^2)"
+
+hist_all.base_geometry = histogram
+hist_all.hist_species_names = beam plasma
+hist_all.hist_num_bins = 200 200
+hist_all.hist_bins_lo = -10 -10
+hist_all.hist_bins_hi = 10 10
+# test all quantities available for the histogram
+hist_all.hist_function = log(abs(x+2*y+z+ion_lev))
+hist_all.hist_function2 = log(abs(3*ux+uy+uz))
+hist_all.hist_weight = log(abs(w+ga_psi))*log(abs(3*ux+uy+uz))*log(abs(x+2*y+z+ion_lev))

src/Hipace.H

@@ -365,7 +365,6 @@ private:
     amrex::Vector< CoulombCollision > m_all_collisions;
 
     void InitDiagnostics (const int step, const amrex::Real time, const bool is_last_step);
-    void FillFieldDiagnostics (const int current_N_level, int islice);
     void FillBeamDiagnostics (const int step, const amrex::Real time, const bool is_last_step);
     void WriteDiagnostics (const int step, const amrex::Real time, const bool is_last_step);
     void FlushDiagnostics ();

src/Hipace.cpp

@@ -813,8 +813,13 @@ Hipace::SolveOneSlice (int islice, int step, bool is_first_step, bool is_last_st
     // get laser insitu diagnostics
     m_multi_laser.InSituComputeDiags(step, islice, m_physical_time, is_last_step);
 
-    // copy fields (and laser) to diagnostic array
-    FillFieldDiagnostics(current_N_level, islice);
+    // copy fields, laser, plasma and beam to diagnostic array
+    m_diags.FillDiagnostics(
+        islice, current_N_level,
+        m_fields, m_multi_laser,
+        m_multi_plasma, m_multi_beam,
+        m_3D_geom
+    );
 
     // plasma field ionization
     for (int lev=0; lev<current_N_level; ++lev) {
@@ -1317,16 +1322,6 @@ Hipace::InitDiagnostics (const int step, const amrex::Real time, const bool is_l
     m_diags.ResizeFDiagFAB(m_3D_geom, m_multi_laser.GetLaserGeom(), step, time, is_last_step);
 }
 
-void
-Hipace::FillFieldDiagnostics (const int current_N_level, int islice)
-{
-    for (auto& fd : m_diags.getFieldData()) {
-        if (fd.m_has_field) {
-            m_fields.Copy(current_N_level, islice, fd, m_3D_geom, m_multi_laser);
-        }
-    }
-}
-
 void
 Hipace::FillBeamDiagnostics (const int step, const amrex::Real time, const bool is_last_step)
 {
@@ -1344,7 +1339,7 @@ Hipace::WriteDiagnostics (const int step, const amrex::Real time, const bool is_
 {
 #ifdef HIPACE_USE_OPENPMD
     if (m_diags.hasAnyFieldOutput(step, time, is_last_step)) {
-        m_openpmd_writer.WriteFieldDiagnostics(m_diags.getFieldData(),
+        m_openpmd_writer.WriteFieldDiagnostics(m_diags.getDiagData(),
             m_multi_laser, m_physical_time, step);
     }
 

src/diagnostics/Diagnostic.H

@@ -15,38 +15,63 @@
 
 #include <vector>
 
+class Fields;
+class MultiBeam;
+class MultiLaser;
+class MultiPlasma;
+
 
 /** \brief This struct holds data for one field diagnostic on one MR level */
-struct FieldDiagnosticData
+struct DiagnosticData
 {
     std::string m_diag_name; /**< name used for input parameters and in the output */
-    enum struct geom_type {
+    enum struct diag_type {
         field,
-        laser
-    } m_base_geom_type = geom_type::field;  /**< which geometry the diagnostics is based on */
-    int m_level = 0; /**< MR level when field_geom is used */
-    int m_slice_dir; /**< Slicing direction */
+        laser,
+        histogram
+    } m_base_diag_type = diag_type::field;  /**< which geometry the diagnostics is based on */
+    /** Number of iterations between consecutive output dumps.
+     * Default is 0, meaning no output */
+    utils::DiagPeriod m_output_period = "0";
+    bool m_has_output = false; /**< if there is output to write */
+    int m_level = 0; /**< MR level */
+    amrex::IntVect m_remove_axis {0, 0, 0}; /**< which axis should be removed */
+    amrex::Vector<std::string> m_axis_labels {}; /**< axis labels in C order */
+    bool m_integrate_along_z = false; /**< resolve or integrate along zeta */
+    amrex::Geometry m_geom_io; /**< Diagnostics geometry of the openPMD output */
+    amrex::Geometry m_realspace_geom; /**< Diagnostics geometry in simulation realspace */
     bool m_include_ghost_cells = false; /**< if ghost cells are included in output */
     bool m_use_custom_size_lo = false; /**< if a user defined diagnostics size should be used (lo)*/
     bool m_use_custom_size_hi = false; /**< if a user defined diagnostics size should be used (hi)*/
     /** 3D array with lower ends of the diagnostics grid */
     amrex::RealVect m_diag_lo {0., 0., 0.};
     /** 3D array with upper ends of the diagnostics grid */
     amrex::RealVect m_diag_hi {0., 0., 0.};
-    amrex::IntVect m_diag_coarsen; /**< xyz coarsening ratio (positive) */
-    int m_nfields; /**< Number of physical fields to write */
+    amrex::IntVect m_diag_coarsen {1, 1, 1}; /**< xyz coarsening ratio (positive) */
+    int m_nfields = 0; /**< Number of physical fields to write */
     amrex::Vector<std::string> m_comps_output; /**< Component names to Write to output file */
     /** Component indexes to Write to output file */
     amrex::Gpu::Buffer<int> m_comps_output_idx;
-    /** Vector over levels, all fields */
-    amrex::FArrayBox m_F;
+    /** All fields */
+    amrex::FArrayBox m_F_real;
     /** FAB for laser */
-    amrex::BaseFab<amrex::GpuComplex<amrex::Real>> m_F_laser;
-    amrex::Geometry m_geom_io; /**< Diagnostics geometry */
-    bool m_has_field; /**< if there is field output to write */
-    /** Number of iterations between consecutive output dumps.
-     * Default is 0, meaning no output */
-    utils::DiagPeriod m_output_period = "0";
+    amrex::BaseFab<amrex::GpuComplex<amrex::Real>> m_F_complex;
+
+    // histogram
+
+    amrex::Vector<std::string> m_hist_species_names; /**< beam and plasma names for histogram */
+    int m_hist_num_dims = 1; /**< 1D or 2D histogram */
+    amrex::Vector<int> m_hist_num_bins; /**< number of bins per hist axis */
+    amrex::Vector<amrex::Real> m_hist_bins_lo; /**< bounds for each hist axis */
+    amrex::Vector<amrex::Real> m_hist_bins_hi; /**< bounds for each hist axis */
+    amrex::FArrayBox m_hist_gpu_fab; /**< hist slice data cache on GPU */
+    /** parser funcitons for histogram */
+    amrex::Parser m_hist_parser_q1;
+    amrex::ParserExecutor<9> m_hist_exe_q1;
+    amrex::Parser m_hist_parser_q2;
+    amrex::ParserExecutor<9> m_hist_exe_q2;
+    amrex::Parser m_hist_parser_w;
+    amrex::ParserExecutor<9> m_hist_exe_w;
 };
 
 
@@ -66,7 +91,7 @@ public:
     amrex::Vector<std::string>& getBeamNames () { return m_output_beam_names; }
 
     /** \brief return data for all field diagnostics */
-    amrex::Vector<FieldDiagnosticData>& getFieldData () { return m_field_data; }
+    amrex::Vector<DiagnosticData>& getDiagData () { return m_diag_data; }
 
     /** \brief determines if a single field diagnostic has any output on in this time step
      *
@@ -75,7 +100,7 @@ public:
      * \param[in] output_time physical time of the current step
      * \param[in] is_last_step if this is the last time step of the simulation
      */
-    static bool hasFieldOutput (const FieldDiagnosticData& fd, int output_step,
+    static bool hasFieldOutput (const DiagnosticData& fd, int output_step,
                                 amrex::Real output_time, bool is_last_step)
     {
         return (fd.m_nfields > 0) &&
@@ -91,7 +116,7 @@ public:
     [[nodiscard]] bool
     hasAnyFieldOutput (int output_step, amrex::Real output_time, bool is_last_step) const
     {
-        for (const auto& fd : m_field_data) {
+        for (const auto& fd : m_diag_data) {
             if (hasFieldOutput(fd, output_step, output_time, is_last_step)) return true;
         }
         return false;
@@ -138,14 +163,6 @@ public:
     [[nodiscard]] bool
     needsTempIndividual () const;
 
-    /** \brief calculate box which possibly was trimmed in case of slice IO
-     *
-     * \param[in] slice_dir slicing direction
-     * \param[in,out] domain_3d domain box to be possibly trimmed to a slice box
-     * \param[in,out] rbox_3d real box to be possibly trimmed to a slice box
-     */
-    void TrimIOBox (int slice_dir, amrex::Box& domain_3d, amrex::RealBox& rbox_3d);
-
     /** \brief resizes the FArrayBox of the diagnostics to the currently calculated box
      *
      * \param[in] field_geom geometry of the full simulation domain for all field MR levels
@@ -158,13 +175,38 @@ public:
                          amrex::Geometry const& laser_geom, int output_step,
                          amrex::Real output_time, bool is_last_step);
 
+    /** \brief calculate where data from slice islice will go to using zeroth-order interpolation
+     * returns if this slice is needed and the destination slice index
+     *
+     * \param[in] fd diag data object
+     * \param[in] islice current slice
+     * \param[in] geom3d 3d domain used for the z direction
+     */
+    static std::pair<bool, int>
+    ReverseShapeFactor (const DiagnosticData& fd, int islice, const amrex::Geometry& geom3d);
+
+    /** \brief Fill diagnostic arrays with the simulation data available on the current slice
+     *
+     * \param[in] islice current slice
+     * \param[in] current_N_level number of MR levels active on the current slice
+     * \param[in] fields Fields object
+     * \param[in] lasers MultiLaser object
+     * \param[in] plasmas MultiPlasma object
+     * \param[in] beams MultiBeam object
+     * \param[in] field_geom geometry for all field MR levels
+     */
+    void
+    FillDiagnostics (int islice, int current_N_level,
+                     Fields& fields, MultiLaser& lasers,
+                     MultiPlasma& plasmas, MultiBeam& beams,
+                     const amrex::Vector<amrex::Geometry>& field_geom);
+
 private:
     amrex::Vector<std::string> m_output_beam_names; /**< Component names to Write to output file */
     /** Number of iterations between consecutive output dumps.
      * Default is 0, meaning no output */
     utils::DiagPeriod m_beam_output_period = "0";
-    amrex::Vector<FieldDiagnosticData> m_field_data; /**< individual field diag data */
-    bool m_initialized = false; /**< if this object is fully initialized */
+    amrex::Vector<DiagnosticData> m_diag_data; /**< individual field diag data */
 };
 
 #endif // DIAGNOSTIC_H_