Add beam eigenemittances to reduced diagnostics.

docs/source/dataanalysis/dataanalysis.rst

@@ -84,7 +84,7 @@ The code writes out the values in an ASCII file prefixed ``reduced_beam_characte
 * ``sig_px``, ``sig_py``, ``sig_pt``
     Standard deviation of the particle momentum deviations (energy difference for ``pt``) normalized by the magnitude of the reference particle momentum (unit: dimensionless)
 * ``emittance_x``, ``emittance_y``, ``emittance_t``
-    Normalized rms beam emittance (unit: meter)
+    Unnormalized rms beam emittances (unit: meter)
 * ``alpha_x``, ``alpha_y``, ``alpha_t``
     Courant-Snyder (Twiss) alpha (unit: dimensionless).  Transverse Twiss functions are calculated after removing correlations with particle energy.
 * ``beta_x``, ``beta_y``, ``beta_t``
@@ -93,8 +93,11 @@ The code writes out the values in an ASCII file prefixed ``reduced_beam_characte
     Horizontal and vertical dispersion (unit: meter)
 * ``dispersion_px``, ``dispersion_py``
     Derivative of horizontal and vertical dispersion (unit: dimensionless)
+* ``emittance_xn``, ``emittance_yn``, ``emittance_tn``
+    Normalized rms beam emittances (unit: meter)
 * ``emittance_1``, ``emittance_2``, ``emittance_3``
     Normalized rms beam eigenemittances (aka mode emittances) (unit: meter)
+    These three diagnostics are written optionally if the flag eigenemittances = True.
 * ``charge``
     Total beam charge (unit: Coulomb)
 

docs/source/usage/parameters.rst

@@ -683,8 +683,8 @@ Diagnostics and output
   This option is ignored for the openPMD output elements (remove them from the lattice to disable).
 
 * ``diag.slice_step_diagnostics`` (``boolean``, optional, default: ``false``)
-  By default, diagnostics is performed at the beginning and end of the simulation.
-  Enabling this flag will write diagnostics every step and slice step
+  By default, diagnostics are computed and written at the beginning and end of the simulation.
+  Enabling this flag will write diagnostics at every step and slice step.
 
 * ``diag.file_min_digits`` (``integer``, optional, default: ``6``)
     The minimum number of digits used for the step number appended to the diagnostic file names.
@@ -694,6 +694,10 @@ Diagnostics and output
   Diagnostics for particles lost in apertures, stored as ``diags/openPMD/particles_lost.*`` at the end of the simulation.
   See the ``beam_monitor`` element for backend values.
 
+* ``diag.eigenemittances`` (``boolean``, optional, default: ``false``)
+  If this flag is enabled, the 3 eigenemittances of the 6D beam distribution are computed and written as diagnostics.
+  This flag is disabled by default to reduce computational cost.
+
 
 .. _running-cpp-parameters-diagnostics-insitu:
 

docs/source/usage/python.rst

@@ -156,6 +156,13 @@ Collective Effects & Overall Simulation Parameters
       Diagnostics for particles lost in apertures.
       See the ``BeamMonitor`` element for backend values.
 
+   .. py:property:: eigenemittances
+
+      Enable (``True``) or disable (``False``) output of eigenemittances at every slice step in elements  (default: ``False``).
+
+      If this flag is enabled, the 3 eigenemittances of the 6D beam distribution are computed and written as diagnostics.
+      This flag is disabled by default to reduce computational cost.
+
    .. py:method:: init_grids()
 
       Initialize AMReX blocks/grids for domain decomposition & space charge mesh.

src/particles/diagnostics/CovarianceMatrixMath.H

@@ -11,6 +11,7 @@
 #define COVARIANCE_MATRIX_MATH_H
 
 #include <ablastr/constant.H>
+#include <ablastr/warn_manager/WarnManager.H>
 
 #include <AMReX_Array.H>
 #include <AMReX_REAL.H>
@@ -62,6 +63,15 @@ namespace impactx::diagnostics
         amrex::ParticleReal tol = 1.0e-12;  //allow for roundoff error
         if (discriminant > tol) {
 
+            ablastr::warn_manager::WMRecordWarning(
+                    "Impactx::diagnostics::CubicRootsTrig",
+                    "Polynomial appearing in CubicRootsTrig has one or more complex "
+                    "(non-real) roots.  Only the real part is returned.  This "
+                    "suggests a loss of numerical precision in computation of the "
+                    "eigenemittances.  Treat eigenemittance values with caution.",
+                    ablastr::warn_manager::WarnPriority::medium
+            );
+
            std::cout << "Polynomial in CubicRoots has one or more complex roots." << "\n";
 
         } else if (Q == 0.0_prt) {  // Special case of a triple root

src/particles/diagnostics/DiagnosticOutput.cpp

@@ -43,7 +43,14 @@ namespace impactx::diagnostics
             if (otype == OutputType::PrintRefParticle) {
                 file_handler << "step s beta gamma beta_gamma x y z t px py pz pt\n";
             } else if (otype == OutputType::PrintReducedBeamCharacteristics) {
-                file_handler << "step" << " " << "s" << " "
+
+                // determine whether to output eigenemittances
+                amrex::ParmParse pp_diag("diag");
+                bool compute_eigenemittances = false;
+                pp_diag.queryAdd("eigenemittances", compute_eigenemittances);
+
+                if (compute_eigenemittances) {
+                   file_handler << "step" << " " << "s" << " "
                              << "x_mean" << " " << "x_min" << " " << "x_max" << " "
                              << "y_mean" << " " << "y_min" << " " << "y_max" << " "
                              << "t_mean" << " " << "t_min" << " " << "t_max" << " "
@@ -57,9 +64,29 @@ namespace impactx::diagnostics
                              << "beta_x" << " " << "beta_y" << " " << "beta_t" << " "
                              << "dispersion_x" << " " << "dispersion_px" << " "
                              << "dispersion_y" << " " << "dispersion_py" << " "
+                             << "emittance_xn" << " " << "emittance_yn" << " " << "emittance_tn" << " "
                              << "emittance_1" << " " << "emittance_2" << " " << "emittance_3" << " "
                              << "charge_C" << " "
                              << "\n";
+                } else {
+                   file_handler << "step" << " " << "s" << " "
+                             << "x_mean" << " " << "x_min" << " " << "x_max" << " "
+                             << "y_mean" << " " << "y_min" << " " << "y_max" << " "
+                             << "t_mean" << " " << "t_min" << " " << "t_max" << " "
+                             << "sig_x" << " " << "sig_y" << " " << "sig_t" << " "
+                             << "px_mean" << " " << "px_min" << " " << "px_max" << " "
+                             << "py_mean" << " " << "py_min" << " " << "py_max" << " "
+                             << "pt_mean" << " " << "pt_min" << " " << "pt_max" << " "
+                             << "sig_px" << " " << "sig_py" << " " << "sig_pt" << " "
+                             << "emittance_x" << " " << "emittance_y" << " " << "emittance_t" << " "
+                             << "alpha_x" << " " << "alpha_y" << " " << "alpha_t" << " "
+                             << "beta_x" << " " << "beta_y" << " " << "beta_t" << " "
+                             << "dispersion_x" << " " << "dispersion_px" << " "
+                             << "dispersion_y" << " " << "dispersion_py" << " "
+                             << "emittance_xn" << " " << "emittance_yn" << " " << "emittance_tn" << " "
+                             << "charge_C" << " "
+                             << "\n";
+                }
             }
         }
 
@@ -93,7 +120,13 @@ namespace impactx::diagnostics
 
             amrex::ParticleReal const s = pc.GetRefParticle().s;
 
-            file_handler << step << " " << s << " "
+            // determine whether to output eigenemittances
+            amrex::ParmParse pp_diag("diag");
+            bool compute_eigenemittances = false;
+            pp_diag.queryAdd("eigenemittances", compute_eigenemittances);
+
+            if (compute_eigenemittances) {
+               file_handler << step << " " << s << " "
                          << rbc.at("x_mean") << " " << rbc.at("x_min") << " " << rbc.at("x_max") << " "
                          << rbc.at("y_mean") << " " << rbc.at("y_min") << " " << rbc.at("y_max") << " "
                          << rbc.at("t_mean") << " " << rbc.at("t_min") << " " << rbc.at("t_max") << " "
@@ -107,8 +140,27 @@ namespace impactx::diagnostics
                          << rbc.at("beta_x") << " " << rbc.at("beta_y") << " " << rbc.at("beta_t") << " "
                          << rbc.at("dispersion_x") << " " << rbc.at("dispersion_px") << " "
                          << rbc.at("dispersion_y") << " " << rbc.at("dispersion_py") << " "
+                         << rbc.at("emittance_xn") << " " << rbc.at("emittance_yn") << " " << rbc.at("emittance_tn") << " "
                          << rbc.at("emittance_1") << " " << rbc.at("emittance_2") << " " << rbc.at("emittance_3") << " "
                          << rbc.at("charge_C") << "\n";
+            } else {
+               file_handler << step << " " << s << " "
+                         << rbc.at("x_mean") << " " << rbc.at("x_min") << " " << rbc.at("x_max") << " "
+                         << rbc.at("y_mean") << " " << rbc.at("y_min") << " " << rbc.at("y_max") << " "
+                         << rbc.at("t_mean") << " " << rbc.at("t_min") << " " << rbc.at("t_max") << " "
+                         << rbc.at("sig_x") << " " << rbc.at("sig_y") << " " << rbc.at("sig_t") << " "
+                         << rbc.at("px_mean") << " " << rbc.at("px_min") << " " << rbc.at("px_max") << " "
+                         << rbc.at("py_mean") << " " << rbc.at("py_min") << " " << rbc.at("py_max") << " "
+                         << rbc.at("pt_mean") << " " << rbc.at("pt_min") << " " << rbc.at("pt_max") << " "
+                         << rbc.at("sig_px") << " " << rbc.at("sig_py") << " " << rbc.at("sig_pt") << " "
+                         << rbc.at("emittance_x") << " " << rbc.at("emittance_y") << " " << rbc.at("emittance_t") << " "
+                         << rbc.at("alpha_x") << " " << rbc.at("alpha_y") << " " << rbc.at("alpha_t") << " "
+                         << rbc.at("beta_x") << " " << rbc.at("beta_y") << " " << rbc.at("beta_t") << " "
+                         << rbc.at("dispersion_x") << " " << rbc.at("dispersion_px") << " "
+                         << rbc.at("dispersion_y") << " " << rbc.at("dispersion_py") << " "
+                         << rbc.at("emittance_xn") << " " << rbc.at("emittance_yn") << " " << rbc.at("emittance_tn") << " "
+                         << rbc.at("charge_C") << "\n";
+            }
         } // if( otype == OutputType::PrintReducedBeamCharacteristics)
 
         // TODO: add as an option to the monitor element

src/particles/diagnostics/EmittanceInvariants.cpp

@@ -59,7 +59,7 @@ namespace impactx::diagnostics
         // a change of sign.
         for (int i = 1; i < 7; i++) {
             for (int j = 1; j < 7; j++) {
-                if (j % 2 == 1) {
+                if (j % 2 != 0) {
                    S1(i,j) = -Sigma(i,j+1); // if j is odd
                 }
                 else {

src/particles/diagnostics/ReducedBeamCharacteristics.cpp

@@ -305,19 +305,22 @@ namespace impactx::diagnostics
         amrex::ParticleReal const alpha_y = - ypy_d / emittance_yd;
         amrex::ParticleReal const alpha_t = - tpt / emittance_t;
 
-        // Calculate eigenemittances (optional)
-        amrex::Array2D<amrex::ParticleReal, 1, 6, 1, 6> Sigma;
-        amrex::ParticleReal emittance_1 = emittance_x * bg;
-        amrex::ParticleReal emittance_2 = emittance_y * bg;
-        amrex::ParticleReal emittance_3 = emittance_t * bg;
+        // Calculate normalized emittances
+        amrex::ParticleReal emittance_xn = emittance_x * bg;
+        amrex::ParticleReal emittance_yn = emittance_y * bg;
+        amrex::ParticleReal emittance_tn = emittance_t * bg;
 
-        // Parse the diagnostic parameters
+        // Determine whether to calculate eigenemittances, and initialize
         amrex::ParmParse pp_diag("diag");
         bool compute_eigenemittances = false;
         pp_diag.queryAdd("eigenemittances", compute_eigenemittances);
+        amrex::ParticleReal emittance_1 = emittance_xn;
+        amrex::ParticleReal emittance_2 = emittance_yn;
+        amrex::ParticleReal emittance_3 = emittance_tn;
 
         if (compute_eigenemittances) {
            // Store the covariance matrix in dynamical variables:
+           amrex::Array2D<amrex::ParticleReal, 1, 6, 1, 6> Sigma;
            Sigma(1,1) = x_ms;
            Sigma(1,2) = xpx * bg;
            Sigma(1,3) = xy;
@@ -400,9 +403,14 @@ namespace impactx::diagnostics
         data["dispersion_y"] = dispersion_y;
         data["dispersion_py"] = dispersion_py;
         data["charge_C"] = charge;
-        data["emittance_1"] = emittance_1;
-        data["emittance_2"] = emittance_2;
-        data["emittance_3"] = emittance_3;
+        data["emittance_xn"] = emittance_xn;
+        data["emittance_yn"] = emittance_yn;
+        data["emittance_tn"] = emittance_tn;
+        if (compute_eigenemittances) {
+           data["emittance_1"] = emittance_1;
+           data["emittance_2"] = emittance_2;
+           data["emittance_3"] = emittance_3;
+        }
 
         return data;
     }

tests/python/test_transformation.py

@@ -81,7 +81,7 @@ def test_transformation():
     for key, val in rbc_s0.items():
         if not np.isclose(val, rbc_s[key], rtol=rtol, atol=atol):
             print(f"initial[{key}]={val}, final[{key}]={rbc_s[key]} not equal")
-            assert np.isclose(val, rbc_s[key], rtol=rtol, atol=atol)
+        assert np.isclose(val, rbc_s[key], rtol=rtol, atol=atol)
     # assert that the t-based beam is different, at least in the following keys:
     large_st_diff_keys = [
         "beta_x",