Plasma Temperature Diagnostics

docs/source/run/parameters.rst

@@ -1018,21 +1018,16 @@ Field diagnostics
     Whether the field diagnostics should include ghost cells.
 
 * ``<diag name> or diagnostic.field_data`` (`string`) optional (default `all`)
-    Names of the fields written to file, separated by a space. The field names need to be ``all``,
-    ``none`` or a subset of ``ExmBy EypBx Ez Bx By Bz Psi``. For the predictor-corrector solver,
-    additionally ``jx jy jz rhomjz`` are available, which are the current and charge densities of the
-    plasma and the beam, with ``rhomjz`` equal to :math:`\rho-j_z/c`.
-    For the explicit solver, the current and charge densities of the beam and
-    for all plasmas are separated: ``jx_beam jy_beam jz_beam`` and ``jx jy rhomjz`` are available.
-    If ``rho`` is explicitly mentioned as ``field_data``, it is deposited by the plasma
-    to be available as a diagnostic. Similarly if ``rho_<plasma name>`` is explicitly mentioned,
-    the charge density of that plasma species will be separately available as a diagnostic.
-    When a laser pulse is used, the laser complex envelope ``laserEnvelope`` is available
-    in the ``laser`` base geometry.
-    The plasma proper density (n/gamma) is then also accessible via ``chi``.
-    A field can be removed from the list, for example, after it has been included through ``all``,
-    by adding ``remove_<field name>`` after it has been added. If a field is added and removed
-    multiple times, the last occurrence takes precedence.
+    Specifies the fields to be written to file, separated by a space. The field names can be:
+
+    * ``all``: Includes all available fields.
+    * ``none``: Excludes all fields.
+    * A subset of the following: ``ExmBy``, ``EypBx``, ``Ez``, ``Bx``, ``By``, ``Bz``, ``Psi``.
+    * Specific to the Predictor-Corrector solver: ``jx``, ``jy``, ``jz``, and ``rhomjz``, which correspond to the current and charge densities of the plasma and beam (``rhomjz`` is defined as :math:`\rho-j_z/c`).
+    * Specific to the Explicit solver: separate current and charge densities for the beam (``jx_beam``, ``jy_beam``, ``jz_beam``) and plasma (``jx``, ``jy``, and ``rhomjz``).
+    * Plasma diagnostics: ``rho`` (total charge density) is always available. Per-species diagnostics are also available: ``rho_<plasma name>`` (charge density of the species); ``w_<plasma name>`` (particle weights of the species); and momentum components ``ux_<plasma name>``, ``uy_<plasma name>``, ``uz_<plasma name>``, ``ux^2_<plasma name>``, etc.
+    * Laser diagnostics, when a laser pulse is used: ``laserEnvelope`` (the complex envelope of the laser in the ``laser`` base geometry) and ``chi`` (plasma proper density :math:`n/\gamma`).
+    * Fields can be added or removed from the list dynamically: to remove a field after including ``all``, use ``remove_<field name>``. If a field is added and removed multiple times, the last occurrence takes precedence.
 
 * ``<diag name> or diagnostic.patch_lo`` (3 `float`) optional (default `-infinity -infinity -infinity`)
     Lower limit for the diagnostic grid.
@@ -1046,9 +1041,14 @@ Field diagnostics
     If ``rho`` is explicitly mentioned in ``diagnostic.field_data``, then the default will become `1`.
 
 * ``hipace.deposit_rho_individual`` (`bool`) optional (default `0`)
-    This option works similar to ``hipace.deposit_rho``,
-    however the charge density from every plasma species will be deposited into individual fields
-    that are accessible as ``rho_<plasma name>`` in ``diagnostic.field_data``.
+    This option works similarly to ``hipace.deposit_rho``,
+    but the charge density from every plasma species will be deposited into individual fields
+    accessible as ``rho_<plasma name>`` in ``diagnostic.field_data``.
+
+* ``hipace.deposit_temp_individual`` (`bool`) optional (default `0`)
+    The weights, momentum, and their squares from every plasma species
+    will be deposited into individual fields accessible as ``w``, ``ux_<plasma name>`` or
+    ``ux^2_<plasma name>`` (similarly for ``uy`` and ``uz``) in ``diagnostic.field_data``.
 
 In-situ diagnostics
 ^^^^^^^^^^^^^^^^^^^

examples/laser_ionization/analysis_laser_ionization.py

@@ -9,7 +9,6 @@
 import numpy as np
 import math
 from openpmd_viewer import OpenPMDTimeSeries
-import statistics
 import scipy.constants as scc
 
 import sys
@@ -48,24 +47,26 @@
 n0 = nc / 10000
 
 iteration = 0
+izmax = 10
 
 rho_elec_linear, _ = ts_linear.get_field(field='rho_elec', coord='z', iteration=iteration)
 rho_elec_mean_linear = np.mean(rho_elec_linear, axis=(1, 2))
-rho_average_linear = statistics.mean(rho_elec_mean_linear[0:10])
+rho_average_linear = np.mean(rho_elec_mean_linear[0:izmax])
 fraction_linear = -rho_average_linear / scc.e / n0
 
 rho_elec_circular, _ = ts_circular.get_field(field='rho_elec', coord='z', iteration=iteration)
 rho_elec_mean_circular = np.mean(rho_elec_circular, axis=(1, 2))
-rho_average_circular = statistics.mean(rho_elec_mean_circular[0:10]) #average over a thickness in the ionized region
+rho_average_circular = np.mean(rho_elec_mean_circular[0:izmax]) #average over a thickness in the ionized region
 fraction_circular = -rho_average_circular / scc.e / n0
 
 fraction_warpx_linear = 0.41014984 # result from WarpX simulation
 fraction_warpx_circular = 0.502250841 # result from WarpX simulation
 
-relative_diff_linear = np.abs( ( fraction_linear - fraction_warpx_linear ) / fraction_warpx_linear )
-relative_diff_circular = np.abs( ( fraction_circular - fraction_warpx_circular ) / fraction_warpx_circular )
+error_fraction_linear = np.abs( ( fraction_linear - fraction_warpx_linear ) / fraction_warpx_linear )
+error_fraction_circular = np.abs( ( fraction_circular - fraction_warpx_circular ) / fraction_warpx_circular )
 
-tolerance = 0.15
+tolerance_higher = 0.15
+tolerance_lower = 0.001
 print(f"fraction_warpx_linear = {fraction_warpx_linear}")
 print(f"fraction_hipace_linear = {fraction_linear}")
 print(f"fraction_warpx_circular = {fraction_warpx_circular}")
@@ -74,32 +75,62 @@
 # in-situ diagnostics for calculation of the temperature in all directions in circular polarization
 insitu_path_linear = f'./{args.third}/reduced_elec.0000.txt'
 all_data_linear = diag.read_file(insitu_path_linear)
-Tx2_l = all_data_linear['[ux^2]'][0,0]*scc.m_e*scc.c**2/scc.e
-Ty2_l = all_data_linear['[uy^2]'][0,0]*scc.m_e*scc.c**2/scc.e
-Tz2_l = all_data_linear['[uz^2]'][0,0]*scc.m_e*scc.c**2/scc.e
-temp_eV_linear = 1./3*(Tx2_l+Ty2_l+Tz2_l)
+Tx2 = all_data_linear['[ux^2]'][0,0]*scc.m_e*scc.c**2/scc.e
+Ty2 = all_data_linear['[uy^2]'][0,0]*scc.m_e*scc.c**2/scc.e
+Tz2 = all_data_linear['[uz^2]'][0,0]*scc.m_e*scc.c**2/scc.e
+temp_eV_linear = 1./3*(Tx2+Ty2+Tz2)
 
 insitu_path_circular = f'./{args.fourth}/reduced_elec.0000.txt'
 all_data_circular = diag.read_file(insitu_path_circular)
-Tx2_c = all_data_circular['[ux^2]'][0,0]*scc.m_e*scc.c**2/scc.e
-Ty2_c = all_data_circular['[uy^2]'][0,0]*scc.m_e*scc.c**2/scc.e
-Tz2_c = all_data_circular['[uz^2]'][0,0]*scc.m_e*scc.c**2/scc.e
-temp_eV_circular = 1./3*(Tx2_c+Ty2_c+Tz2_c)
+Tx2 = all_data_circular['[ux^2]'][0,0]*scc.m_e*scc.c**2/scc.e
+Ty2 = all_data_circular['[uy^2]'][0,0]*scc.m_e*scc.c**2/scc.e
+Tz2 = all_data_circular['[uz^2]'][0,0]*scc.m_e*scc.c**2/scc.e
+temp_eV_circular = 1./3*(Tx2+Ty2+Tz2)
+
+# calculation of temperature with OpenPMD-viewer diagnostics
+wf, _ = ts_linear.get_field(field='w_elec', iteration=iteration)
+sum_wf = np.sum(wf, axis=(1,2))[0:izmax]
+ux2_linear, _ = ts_linear.get_field(field='ux^2_elec', iteration=iteration)
+uy2_linear, _ = ts_linear.get_field(field='uy^2_elec', iteration=iteration)
+uz2_linear, _ = ts_linear.get_field(field='uz^2_elec', iteration=iteration)
+ux2_average_linear = np.mean(np.sum(ux2_linear * wf, axis=(1, 2))[0:izmax] / sum_wf)
+uy2_average_linear = np.mean(np.sum(uy2_linear * wf, axis=(1, 2))[0:izmax] / sum_wf)
+uz2_average_linear = np.mean(np.sum(uz2_linear * wf, axis=(1, 2))[0:izmax] / sum_wf)
+temp_diags_linear = 1./3*(ux2_average_linear + uy2_average_linear + uz2_average_linear)*scc.m_e*scc.c**2/scc.e
+
+wf, _ = ts_circular.get_field(field='w_elec', iteration=iteration)
+sum_wf = np.sum(wf, axis=(1,2))[0:izmax]
+ux2_circular, _ = ts_circular.get_field(field='ux^2_elec', iteration=iteration)
+uy2_circular, _ = ts_circular.get_field(field='uy^2_elec', iteration=iteration)
+uz2_circular, _ = ts_circular.get_field(field='uz^2_elec', iteration=iteration)
+ux2_average_circular = np.mean(np.sum(ux2_circular * wf, axis=(1, 2))[0:izmax] / sum_wf)
+uy2_average_circular = np.mean(np.sum(uy2_circular * wf, axis=(1, 2))[0:izmax] / sum_wf)
+uz2_average_circular = np.mean(np.sum(uz2_circular * wf, axis=(1, 2))[0:izmax] / sum_wf)
+temp_diags_circular = 1./3*(ux2_average_circular + uy2_average_circular + uz2_average_circular)*scc.m_e*scc.c**2/scc.e
 
 temp_eV_warpx_linear = 1.00286009
 temp_eV_warpx_circular = 9.68224535
 
-print(f"temp_eV_warpx_linear = {temp_eV_warpx_linear}")
-print(f"temp_eV_hipace_linear = {temp_eV_linear}")
+print(f"temperature (eV) WarpX, linear = {temp_eV_warpx_linear}")
+print(f"temperature (eV) HiPACE++, linear with insitu diagnostics = {temp_eV_linear}")
+print(f"temperature (eV) HiPACE++, linear with Open-PMD diagnostics = {temp_diags_linear}")
 
-print(f"temp_eV_warpx_circular = {temp_eV_warpx_circular}")
-print(f"temp_eV_hipace_circular = {temp_eV_circular}")
+print(f"temperature (eV) WarpX, circular = {temp_eV_warpx_circular}")
+print(f"temperature (eV) HiPACE++, circular with insitu diagnostics = {temp_eV_circular}")
+print(f"temperature (eV) HiPACE++, circular with Open-PMD diagnostics = {temp_diags_circular}")
 
-relative_diff_temp_linear = np.abs( ( temp_eV_linear - temp_eV_warpx_linear ) / temp_eV_warpx_linear )
-relative_diff_temp_circular = np.abs( ( temp_eV_circular - temp_eV_warpx_circular ) / temp_eV_warpx_circular )
+# Error of temperature calculation between insitu-diagnostics and Open-PMD diagnostics in HiPACE++
+error_h_h_linear = np.abs( ( temp_eV_linear - temp_diags_linear ) / temp_eV_linear )
+error_h_h_circular = np.abs( ( temp_eV_circular - temp_diags_circular ) / temp_eV_circular )
 
-assert ( (relative_diff_linear < tolerance) and \
-         (relative_diff_circular < tolerance) and \
-         (relative_diff_temp_linear < tolerance) and \
-         (relative_diff_temp_circular < tolerance)), \
+# Error of temperature calculation between HiPACE++ and WarpX
+error_h_w_linear = np.abs( ( temp_eV_linear - temp_eV_warpx_linear ) / temp_eV_warpx_linear )
+error_h_w_circular = np.abs( ( temp_eV_circular - temp_eV_warpx_circular ) / temp_eV_warpx_circular )
+
+assert ( (error_fraction_linear < tolerance_higher) and \
+         (error_fraction_circular < tolerance_higher) and \
+         (error_h_h_linear < tolerance_lower) and \
+         (error_h_h_linear < tolerance_lower) and \
+         (error_h_w_linear < tolerance_higher) and \
+         (error_h_w_circular < tolerance_higher)), \
          'Test laser_ionization did not pass'

examples/laser_ionization/inputs_laser_ionization

@@ -4,29 +4,26 @@ hipace.verbose = 3
 
 amr.n_cell = 50 50 50
 
-
 my_constants.kp_inv = 10.e-6
 
 my_constants.nc = 1.75e27
 my_constants.n0 = nc/10000
 
-my_constants.a0 = 0.00885126
-
 my_constants.w0_um = 1.e8
 my_constants.tau_fs = 30/1.17741 
 my_constants.lambda0_nm = 800
 
 hipace.file_prefix = laser_ionization.1Rank
 hipace.do_tiling = 0
 hipace.deposit_rho_individual = 1
+hipace.deposit_temp_individual = 1
 
 geometry.prob_lo     = -10.*kp_inv   -10.*kp_inv   -15.*kp_inv #box shifted
 geometry.prob_hi     =  10.*kp_inv    10.*kp_inv    5.*kp_inv
 
 lasers.names = laser
 lasers.lambda0 = lambda0_nm * 1.e-9
 
-
 laser.a0 = a0
 laser.position_mean = 0. 0. 0
 laser.w0 = w0_um*1.e-6

src/Hipace.H

@@ -229,6 +229,8 @@ public:
     inline static bool m_deposit_rho = false;
     /** Whether to deposit rho for every individual plasma for diagnostics */
     inline static bool m_deposit_rho_individual = false;
+    /** Whether to deposit temperature for every individual plasma for diagnostics */
+    inline static bool m_deposit_temp_individual = false;
     /** Whether to interpolate the neutralizing background to MR levels 1 and 2 instead of depositing it */
     inline static bool m_interpolate_neutralizing_background = false;
     /** Whether to use tiling for particle operations */

src/Hipace.cpp

@@ -117,6 +117,8 @@ Hipace::Hipace () :
     queryWithParser(pph, "deposit_rho", m_deposit_rho);
     m_deposit_rho_individual = m_diags.needsRhoIndividual();
     queryWithParser(pph, "deposit_rho_individual", m_deposit_rho_individual);
+    m_deposit_temp_individual = m_diags.needsTempIndividual();
+    queryWithParser(pph, "deposit_temp_individual", m_deposit_temp_individual);
     queryWithParser(pph, "interpolate_neutralizing_background",
         m_interpolate_neutralizing_background);
     bool do_mfi_sync = false;
@@ -595,8 +597,6 @@ Hipace::SolveOneSlice (int islice, int step)
         m_multi_beam.ReorderParticles( WhichBeamSlice::This, step, m_slice_geom[0]);
     }
 
-    m_multi_plasma.InSituComputeDiags(step, islice, m_max_step, m_physical_time, m_max_time);
-
     if (m_N_level > 1) {
         m_multi_beam.TagByLevel(current_N_level, m_3D_geom, WhichSlice::This);
         m_multi_plasma.TagByLevel(current_N_level, m_3D_geom);
@@ -613,6 +613,15 @@ Hipace::SolveOneSlice (int islice, int step)
     // write laser aabs into fields MultiFab
     m_multi_laser.UpdateLaserAabs(islice, current_N_level, m_fields, m_3D_geom);
 
+    // has to be after aabs writing
+    m_multi_plasma.InSituComputeDiags(step, islice, m_max_step, m_physical_time, m_max_time);
+
+    // deposit temperature
+    for (int lev=0; lev<current_N_level; ++lev) {
+        // deposit w, ux, uy, uz, ux2, uy2 and uz2 for all plasmas
+        m_multi_plasma.DoDepositTemperature(m_fields, m_3D_geom, lev);
+    }
+
     // deposit current
     for (int lev=0; lev<current_N_level; ++lev) {
         if (m_explicit) {

src/diagnostics/Diagnostic.H

@@ -138,6 +138,10 @@ public:
      */
     bool needsRhoIndividual () const;
 
+    /** \brief determines if the temperature parameters for every individual plasma are requested as a diagnostic
+     */
+    bool needsTempIndividual () const;
+
     /** \brief calculate box which possibly was trimmed in case of slice IO
      *
      * \param[in] slice_dir slicing direction

src/diagnostics/Diagnostic.cpp

@@ -120,6 +120,25 @@ Diagnostic::needsRhoIndividual () const {
     return false;
 }
 
+bool
+Diagnostic::needsTempIndividual () const {
+    amrex::ParmParse ppd("diagnostic");
+    for (auto& fd : m_field_data) {
+        amrex::ParmParse pp(fd.m_diag_name);
+        amrex::Vector<std::string> comps{};
+        queryWithParserAlt(pp, "field_data", comps, ppd);
+        for (auto& c : comps) {
+            // we don't know the names of all the plasmas here so just look for "ux_..."
+            if (c.find("w_") == 0 ||
+                c.find("ux_") == 0 || c.find("uy_") == 0 || c.find("uz_") == 0 ||
+                c.find("ux^2_") == 0 || c.find("uy^2_") == 0 || c.find("uz^2_") == 0) {
+                return true;
+            }
+        }
+    }
+    return false;
+}
+
 void
 Diagnostic::Initialize (int nlev, bool use_laser) {
     amrex::ParmParse ppd("diagnostic");

src/fields/Fields.cpp

@@ -96,6 +96,12 @@ Fields::AllocData (
                     Comps[isl].multi_emplace(N_Comps, "rho_" + plasma_name);
                 }
             }
+            if (Hipace::m_deposit_temp_individual) {
+                for (auto& plasma_name : Hipace::GetInstance().m_multi_plasma.GetNames()) {
+                    Comps[isl].multi_emplace(N_Comps, "w_" + plasma_name, "ux_" + plasma_name, "uy_" + plasma_name,
+                    "uz_" + plasma_name, "ux^2_" + plasma_name, "uy^2_" + plasma_name, "uz^2_" + plasma_name);
+                }
+            }
             if (Hipace::m_do_beam_jz_minus_rho) {
                 Comps[isl].multi_emplace(N_Comps, "rhomjz_beam");
             }
@@ -144,6 +150,12 @@ Fields::AllocData (
                     Comps[isl].multi_emplace(N_Comps, "rho_" + plasma_name);
                 }
             }
+            if (Hipace::m_deposit_temp_individual) {
+                for (auto& plasma_name : Hipace::GetInstance().m_multi_plasma.GetNames()) {
+                    Comps[isl].multi_emplace(N_Comps, "w_" + plasma_name, "ux_" + plasma_name, "uy_" + plasma_name,
+                    "uz_" + plasma_name, "ux^2_" + plasma_name, "uy^2_" + plasma_name, "uz^2_" + plasma_name);
+                }
+            }
 
             isl = WhichSlice::Previous;
             Comps[isl].multi_emplace(N_Comps, "Bx", "By", "jx", "jy");
@@ -583,6 +595,12 @@ Fields::InitializeSlices (int lev, int islice, const amrex::Vector<amrex::Geomet
             setVal(0., lev, WhichSlice::This, "rho_" + plasma_name);
         }
     }
+    if (Hipace::m_deposit_temp_individual) {
+        for (auto& plasma_name : Hipace::GetInstance().m_multi_plasma.GetNames()) {
+            setVal(0., lev, WhichSlice::This, "w_" + plasma_name, "ux_" + plasma_name, "uy_" + plasma_name,
+            "uz_" + plasma_name, "ux^2_" + plasma_name, "uy^2_" + plasma_name, "uz^2_" + plasma_name);
+        }
+    }
 }
 
 void

src/particles/deposition/CMakeLists.txt

@@ -8,6 +8,7 @@
 target_sources(HiPACE
   PRIVATE
     BeamDepositCurrent.cpp
+    TemperatureDeposition.cpp
     PlasmaDepositCurrent.cpp
     ExplicitDeposition.cpp
 )

src/particles/deposition/TemperatureDeposition.H

@@ -0,0 +1,27 @@
+/* Copyright 2020-2021
+ *
+ * This file is part of HiPACE++.
+ *
+ * Authors: EyaDammak, AlexanderSinn
+ * License: BSD-3-Clause-LBNL
+ */
+#ifndef TEMPERATUREDEPOSITION_H_
+#define TEMPERATUREDEPOSITION_H_
+
+#include "particles/plasma/PlasmaParticleContainer.H"
+#include "fields/Fields.H"
+#include "utils/Constants.H"
+#include "Hipace.H"
+
+/** Depose current of particles in species plasma into the current 2D slice in fields
+ * \param[in] plasma species of which the current is deposited
+ * \param[in,out] fields the general field class, modified by this function
+ * \param[in] gm Geometry of the simulation, to get the cell size etc.
+ * \param[in] lev MR level
+ */
+void
+DepositTemperature (PlasmaParticleContainer& plasma, Fields & fields,
+                amrex::Vector<amrex::Geometry> const& gm, int const lev);
+
+
+#endif //  TEMPERATUREDEPOSITION_H_