Merge pull request #580 from Yin-YinjianZhao/reduced_diags

Add Reduced Diagnostics

18 files changed, 954 insertions, 1 deletions

diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst
index 7db069a6d..536688609 100644
--- a/Docs/source/running_cpp/parameters.rst
+++ b/Docs/source/running_cpp/parameters.rst
@@ -1091,6 +1091,78 @@ Diagnostics and output
     slice diagnostic if there are within the user-defined width from
     the slice region defined by ``slice.dom_lo`` and ``slice.dom_hi``.
 
+* ``warpx.reduced_diags_names`` (`strings`, separated by spaces)
+    The names given by the user of simple reduced diagnostics.
+    Also the names of the output `.txt` files.
+    This reduced diagnostics aims to produce simple outputs
+    of the time history of some physical quantities.
+    If ``warpx.reduced_diags_names`` is not provided in the input file,
+    no reduced diagnostics will be done.
+    This is then used in the rest of the input deck;
+    in this documentation we use `<reduced_diags_name>` as a placeholder.
+
+* ``<reduced_diags_name>.type`` (`string`)
+    The type of reduced diagnostics associated with this `<reduced_diags_name>`.
+    For example, ``ParticleEnergy`` and ``FieldEnergy``.
+    All available types will be descriped below in detail.
+    For all reduced diagnostics,
+    the first and the second columns in the output file are
+    the time step and the corresponding physical time in seconds, respectively.
+
+    * ``ParticleEnergy``
+        This type computes both the total and the mean
+        relativistic particle kinetic energy among all species.
+
+        .. math::
+
+            E_p = \sum_{i=1}^N ( \sqrt{ p_i^2 c^2 + m_0^2 c^4 } - m_0 c^2 ) w_i
+
+        where :math:`p` is the relativistic momentum,
+        :math:`c` is the speed of light,
+        :math:`m_0` is the rest mass,
+        :math:`N` is the number of particles,
+        :math:`w` is the individual particle weight.
+
+        The output columns are
+        total :math:`E_p` of all species,
+        :math:`E_p` of each species,
+        total mean energy :math:`E_p / \sum w_i`,
+        mean enregy of each species.
+
+    * ``FieldEnergy``
+        This type computes the electric and magnetic field energy.
+
+        .. math::
+
+            E_f = \sum [ \varepsilon_0 E^2 / 2 + B^2 / ( 2 \mu_0 ) ] \Delta V
+
+        where
+        :math:`E` is the electric field,
+        :math:`B` is the magnetic field,
+        :math:`\varepsilon_0` is the vacuum permittivity,
+        :math:`\mu_0` is the vacuum permeability,
+        :math:`\Delta V` is the cell volume (or area for 2D),
+        the sum is over all cells.
+
+        The output columns are
+        total field energy :math:`E_f`,
+        :math:`E` field energy,
+        :math:`B` field energy, at mesh refinement levels from 0 to :math:`n`.
+
+* ``<reduced_diags_name>.frequency`` (`int`)
+    The output frequency (every # time steps).
+
+* ``<reduced_diags_name>.path`` (`string`) optional (default `./diags/reducedfiles/`)
+    The path that the output file will be stored.
+
+* ``<reduced_diags_name>.extension`` (`string`) optional (default `txt`)
+    The extension of the output file.
+
+* ``<reduced_diags_name>.separator`` (`string`) optional (default `,`)
+    The separator between row values in the output file.
+    The default separator is comma, i.e. the output file is in
+    the CSV (comma separated value) format.
+
 Lookup tables for QED modules (implementation in progress)
 ----------------------------------------------------------
 Lookup tables store pre-computed values for functions used by the QED modules.
diff --git a/Examples/Tests/reduced_diags/analysis_reduced_diags.py b/Examples/Tests/reduced_diags/analysis_reduced_diags.py
new file mode 100755
index 000000000..f2f8cb673
--- /dev/null
+++ b/Examples/Tests/reduced_diags/analysis_reduced_diags.py
@@ -0,0 +1,88 @@
+#! /usr/bin/env python
+
+# Copyright 2019-2020 Yinjian Zhao
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+# This script tests the reduced diagnostics.
+# The setup is a uniform plasma with electrons and protons.
+# Particle energy and field energy will be outputed
+# using the reduced diagnostics.
+# And they will be compared with the data in the plotfiles.
+
+# Tolerance: 1.0e-8 for particle energy, 1.0e-3 for field energy.
+# The difference of the field energy is relatively large,
+# because fields data in plotfiles are cell-centered,
+# but fields data in reduced diagnostics are staggered.
+
+# Possible running time: ~ 2 s
+
+import sys
+import yt
+import numpy as np
+import scipy.constants as scc
+import csv
+
+fn = sys.argv[1]
+
+ds = yt.load(fn)
+ad = ds.all_data()
+
+# PART1: get results from plotfiles
+
+EPyt = 0.0
+# electron
+px = ad['electrons','particle_momentum_x'].to_ndarray()
+py = ad['electrons','particle_momentum_y'].to_ndarray()
+pz = ad['electrons','particle_momentum_z'].to_ndarray()
+w  = ad['electrons','particle_weight'].to_ndarray()
+EPyt = EPyt + np.sum( (np.sqrt((px**2+py**2+pz**2)*scc.c**2+scc.m_e**2*scc.c**4)-scc.m_e*scc.c**2)*w )
+# proton
+px = ad['protons','particle_momentum_x'].to_ndarray()
+py = ad['protons','particle_momentum_y'].to_ndarray()
+pz = ad['protons','particle_momentum_z'].to_ndarray()
+w  = ad['protons','particle_weight'].to_ndarray()
+EPyt = EPyt + np.sum( (np.sqrt((px**2+py**2+pz**2)*scc.c**2+scc.m_p**2*scc.c**4)-scc.m_p*scc.c**2)*w )
+
+ad = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Ex = ad['Ex'].to_ndarray()
+Ey = ad['Ey'].to_ndarray()
+Ez = ad['Ez'].to_ndarray()
+Bx = ad['Bx'].to_ndarray()
+By = ad['By'].to_ndarray()
+Bz = ad['Bz'].to_ndarray()
+Es = np.sum(Ex**2)+np.sum(Ey**2)+np.sum(Ez**2)
+Bs = np.sum(Bx**2)+np.sum(By**2)+np.sum(Bz**2)
+N  = np.array( ds.domain_width / ds.domain_dimensions )
+dV = N[0]*N[1]*N[2]
+EFyt = 0.5*Es*scc.epsilon_0*dV + 0.5*Bs/scc.mu_0*dV
+
+# PART2: get results from reduced diagnostics
+
+with open('./diags/reducedfiles/EF.txt') as csv_file:
+    csv_reader = csv.reader(csv_file, delimiter=',')
+    line_count = 0
+    for row in csv_reader:
+        if line_count == 2:
+            EF = np.array(row[2]).astype(np.float)
+        line_count += 1
+
+with open('./diags/reducedfiles/EP.txt') as csv_file:
+    csv_reader = csv.reader(csv_file, delimiter=',')
+    line_count = 0
+    for row in csv_reader:
+        if line_count == 2:
+            EP = np.array(row[2]).astype(np.float)
+        line_count += 1
+
+# PART3: print and assert
+
+print('difference of field energy:', abs(EFyt-EF))
+print('tolerance of field energy:', 1.0e-3)
+print('difference of particle energy:', abs(EPyt-EP))
+print('tolerance of particle energy:', 1.0e-10)
+
+assert(abs(EFyt-EF) < 1.0e-3)
+assert(abs(EPyt-EP) < 1.0e-8)
diff --git a/Examples/Tests/reduced_diags/inputs b/Examples/Tests/reduced_diags/inputs
new file mode 100644
index 000000000..75f7a9175
--- /dev/null
+++ b/Examples/Tests/reduced_diags/inputs
@@ -0,0 +1,71 @@
+# 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.coord_sys   =  0            # 0: Cartesian
+geometry.is_periodic =  1    1    1  # Is periodic?
+geometry.prob_lo     = -1.  -1.  -1. # physical domain
+geometry.prob_hi     =  1.   1.   1.
+
+# Algorithms
+algo.current_deposition = esirkepov
+algo.field_gathering = energy-conserving # or momentum-conserving
+warpx.use_filter = 1
+algo.maxwell_fdtd_solver = yee # or ckc
+
+# Interpolation
+# 1: Linear; 2: Quadratic; 3: Cubic.
+interpolation.nox = 1
+interpolation.noy = 1
+interpolation.noz = 1
+
+# CFL
+warpx.cfl = 0.99999
+
+# plot frequency
+amr.plot_int = 200
+
+# Particles
+particles.nspecies = 2
+particles.species_names = electrons protons
+
+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 = gaussian
+protons.ux_th = 0.
+protons.uy_th = 0.
+protons.uz_th = 0.
+
+#################################
+###### REDUCED DIAGS ############
+#################################
+warpx.reduced_diags_names = EP EF
+EP.type = ParticleEnergy
+EP.frequency = 200
+EF.type = FieldEnergy
+EF.frequency = 200
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 2a0923b98..1685c2d2a 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -1082,3 +1082,17 @@ numthreads = 2
 compileTest = 0
 doVis = 0
 compareParticles = 0
+
+[reduced_diags]
+buildDir = .
+inputFile = Examples/Tests/reduced_diags/inputs
+dim = 3
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 2
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/reduced_diags/analysis_reduced_diags.py
diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package
index 710e4399b..12560b49e 100644
--- a/Source/Diagnostics/Make.package
+++ b/Source/Diagnostics/Make.package
@@ -15,5 +15,7 @@ ifeq ($(USE_OPENPMD), TRUE)
   CEXE_sources += WarpXOpenPMD.cpp
 endif
 
+include $(WARPX_HOME)/Source/Diagnostics/ReducedDiags/Make.package
+
 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics
 VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Diagnostics
diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.H b/Source/Diagnostics/ReducedDiags/FieldEnergy.H
new file mode 100644
index 000000000..82fa4b6c4
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.H
@@ -0,0 +1,36 @@
+/* Copyright 2019-2020 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDENERGY_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_FIELDENERGY_H_
+
+#include "ReducedDiags.H"
+#include <fstream>
+
+/**
+ *  This class mainly contains a function that
+ *  computes the field energy.
+ */
+class FieldEnergy : public ReducedDiags
+{
+public:
+
+    /** constructor
+     *  @param[in] rd_name reduced diags names */
+    FieldEnergy(std::string rd_name);
+
+    /** This funciton computes the field energy (EF).
+     *  EF = E eps / 2 + B / mu / 2,
+     *  where E is the electric field,
+     *  B is the magnetic field,
+     *  eps is the vacuum permittivity,
+     *  mu is the vacuum permeability. */
+    virtual void ComputeDiags(int step) override final;
+
+};
+
+#endif
diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp
new file mode 100644
index 000000000..73e6a1c9a
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp
@@ -0,0 +1,139 @@
+/* Copyright 2019-2020 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "FieldEnergy.H"
+#include "WarpX.H"
+#include "WarpXConst.H"
+#include "AMReX_REAL.H"
+#include "AMReX_ParticleReduce.H"
+#include <iostream>
+#include <cmath>
+
+using namespace amrex;
+
+// constructor
+FieldEnergy::FieldEnergy (std::string rd_name)
+: ReducedDiags{rd_name}
+{
+
+    // RZ coordinate is not working
+    #if (defined WARPX_DIM_RZ)
+    AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false,
+        "FieldEnergy reduced diagnostics does not work for RZ coordinate.");
+    #endif
+
+    // get WarpX class object
+    auto & warpx = WarpX::GetInstance();
+
+    // read number of levels
+    int nLevel = 0;
+    ParmParse pp("amr");
+    pp.query("max_level", nLevel);
+    nLevel += 1;
+
+    // resize data array
+    m_data.resize(3*nLevel,0.0);
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        if ( m_IsNotRestart )
+        {
+            // open file
+            std::ofstream ofs;
+            ofs.open(m_path + m_rd_name + "." + m_extension,
+                std::ofstream::out | std::ofstream::app);
+            // write header row
+            ofs << "#";
+            ofs << "[1]step";
+            ofs << m_sep;
+            ofs << "[2]time(s)";
+            for (int lev = 0; lev < nLevel; ++lev)
+            {
+                ofs << m_sep;
+                ofs << "[" + std::to_string(3+3*lev) + "]";
+                ofs << "total(J)lev"+std::to_string(lev);
+                ofs << m_sep;
+                ofs << "[" + std::to_string(4+3*lev) + "]";
+                ofs << "E(J)lev"+std::to_string(lev);
+                ofs << m_sep;
+                ofs << "[" + std::to_string(5+3*lev) + "]";
+                ofs << "B(J)lev"+std::to_string(lev);
+            }
+            ofs << std::endl;
+            // close file
+            ofs.close();
+        }
+    }
+
+}
+// end constructor
+
+// function that computes field energy
+void FieldEnergy::ComputeDiags (int step)
+{
+
+    // Judge if the diags should be done
+    if ( (step+1) % m_freq != 0 ) { return; }
+
+    // get WarpX class object
+    auto & warpx = WarpX::GetInstance();
+
+    // get number of level
+    auto nLevel = warpx.finestLevel() + 1;
+
+    // loop over refinement levels
+    for (int lev = 0; lev < nLevel; ++lev)
+    {
+
+        // get MultiFab data at lev
+        const MultiFab & Ex = warpx.getEfield(lev,0);
+        const MultiFab & Ey = warpx.getEfield(lev,1);
+        const MultiFab & Ez = warpx.getEfield(lev,2);
+        const MultiFab & Bx = warpx.getBfield(lev,0);
+        const MultiFab & By = warpx.getBfield(lev,1);
+        const MultiFab & Bz = warpx.getBfield(lev,2);
+
+        // get cell size
+        Geometry const & geom = warpx.Geom(lev);
+        auto domain_box = geom.Domain();
+        #if (AMREX_SPACEDIM == 2)
+        auto dV = geom.CellSize(0) * geom.CellSize(1);
+        #elif (AMREX_SPACEDIM == 3)
+        auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
+        #endif
+
+        // compute E squared
+        Real tmpx = Ex.norm2(0,geom.periodicity());
+        Real tmpy = Ey.norm2(0,geom.periodicity());
+        Real tmpz = Ez.norm2(0,geom.periodicity());
+        Real Es = tmpx*tmpx + tmpy*tmpy + tmpz*tmpz;
+
+        // compute B squared
+        tmpx = Bx.norm2(0,geom.periodicity());
+        tmpy = By.norm2(0,geom.periodicity());
+        tmpz = Bz.norm2(0,geom.periodicity());
+        Real Bs = tmpx*tmpx + tmpy*tmpy + tmpz*tmpz;
+
+        // save data
+        m_data[lev*3+1] = 0.5 * Es * PhysConst::ep0 * dV;
+        m_data[lev*3+2] = 0.5 * Bs / PhysConst::mu0 * dV;
+        m_data[lev*3+0] = m_data[lev*3+1] + m_data[lev*3+2];
+
+    }
+    // end loop over refinement levels
+
+    /* m_data now contains up-to-date values for:
+     *  [total field energy at level 0,
+     *   electric field energy at level 0,
+     *   magnetic field energy at level 0,
+     *   total field energy at level 1,
+     *   electric field energy at level 1,
+     *   magnetic field energy at level 1,
+     *   ......] */
+
+}
+// end void FieldEnergy::ComputeDiags
diff --git a/Source/Diagnostics/ReducedDiags/Make.package b/Source/Diagnostics/ReducedDiags/Make.package
new file mode 100644
index 000000000..37f76d3d5
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/Make.package
@@ -0,0 +1,14 @@
+CEXE_headers += MultiReducedDiags.H
+CEXE_sources += MultiReducedDiags.cpp
+
+CEXE_headers += ReducedDiags.H
+CEXE_sources += ReducedDiags.cpp
+
+CEXE_headers += ParticleEnergy.H
+CEXE_sources += ParticleEnergy.cpp
+
+CEXE_headers += FieldEnergy.H
+CEXE_sources += FieldEnergy.cpp
+
+INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags
+VPATH_LOCATIONS   += $(WARPX_HOME)/Source/Diagnostics/ReducedDiags
diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H
new file mode 100644
index 000000000..79f487d81
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.H
@@ -0,0 +1,47 @@
+/* Copyright 2019-2020 Maxence Thevenet, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_MULTIREDUCEDDIAGS_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_MULTIREDUCEDDIAGS_H_
+
+#include "ReducedDiags.H"
+#include <vector>
+#include <string>
+#include <memory>
+
+/**
+ * This class holds multiple instances of ReducedDiagnostics, and contains
+ * general functions to initialize, compute, and write these diagnostics
+ * to file.
+ */
+class MultiReducedDiags
+{
+public:
+
+    /// Bool: whether or not reduced diagnostics are activated
+    int m_plot_rd = 0;
+
+    /// names of reduced diagnostics
+    std::vector<std::string> m_rd_names;
+
+    /// m_multi_rd stores a pointer to each reduced diagnostics
+    std::vector<std::unique_ptr<ReducedDiags>> m_multi_rd;
+
+    /// constructor
+    MultiReducedDiags();
+
+    /** Loop over all ReducedDiags and call their ComputeDiags
+     *  @param[in] step current iteration time */
+    void ComputeDiags(int step);
+
+    /** Loop over all ReducedDiags and call their WriteToFile
+     *  @param[in] step current iteration time */
+    void WriteToFile(int step);
+
+};
+
+#endif
diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp
new file mode 100644
index 000000000..71a33924b
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp
@@ -0,0 +1,94 @@
+/* Copyright 2019-2020 Maxence Thevenet, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "ParticleEnergy.H"
+#include "FieldEnergy.H"
+#include "MultiReducedDiags.H"
+#include "AMReX_ParmParse.H"
+#include "AMReX_ParallelDescriptor.H"
+#include <fstream>
+
+using namespace amrex;
+
+// constructor
+MultiReducedDiags::MultiReducedDiags ()
+{
+
+    // read reduced diags names
+    ParmParse pp("warpx");
+    m_plot_rd = pp.queryarr("reduced_diags_names", m_rd_names);
+
+    // if names are not given, reduced diags will not be done
+    if ( m_plot_rd == 0 ) { return; }
+
+    // resize
+    m_multi_rd.resize(m_rd_names.size());
+
+    // loop over all reduced diags
+    for (int i_rd = 0; i_rd < m_rd_names.size(); ++i_rd)
+    {
+
+        ParmParse pp(m_rd_names[i_rd]);
+
+        // read reduced diags type
+        std::string rd_type;
+        pp.query("type", rd_type);
+
+        // match diags
+        if (rd_type.compare("ParticleEnergy") == 0)
+        {
+            m_multi_rd[i_rd].reset
+                ( new ParticleEnergy(m_rd_names[i_rd]));
+        }
+        else if (rd_type.compare("FieldEnergy") == 0)
+        {
+            m_multi_rd[i_rd].reset
+                ( new FieldEnergy(m_rd_names[i_rd]));
+        }
+        else
+        { Abort("No matching reduced diagnostics type found."); }
+        // end if match diags
+
+    }
+    // end loop over all reduced diags
+
+}
+// end constructor
+
+// call functions to compute diags
+void MultiReducedDiags::ComputeDiags (int step)
+{
+    // loop over all reduced diags
+    for (int i_rd = 0; i_rd < m_rd_names.size(); ++i_rd)
+    {
+        m_multi_rd[i_rd] -> ComputeDiags(step);
+    }
+    // end loop over all reduced diags
+}
+// end void MultiReducedDiags::ComputeDiags
+
+// funciton to write data
+void MultiReducedDiags::WriteToFile (int step)
+{
+
+    // Only the I/O rank does
+    if ( !ParallelDescriptor::IOProcessor() ) { return; }
+
+    // loop over all reduced diags
+    for (int i_rd = 0; i_rd < m_rd_names.size(); ++i_rd)
+    {
+
+        // Judge if the diags should be done
+        if ( (step+1) % m_multi_rd[i_rd]->m_freq != 0 ) { return; }
+
+        // call the write to file function
+        m_multi_rd[i_rd]->WriteToFile(step);
+
+    }
+    // end loop over all reduced diags
+}
+// end void MultiReducedDiags::WriteToFile
diff --git a/Source/Diagnostics/ReducedDiags/ParticleEnergy.H b/Source/Diagnostics/ReducedDiags/ParticleEnergy.H
new file mode 100644
index 000000000..d7c60a24e
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/ParticleEnergy.H
@@ -0,0 +1,37 @@
+/* Copyright 2019-2020 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_PARTICLEENERGY_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_PARTICLEENERGY_H_
+
+#include "ReducedDiags.H"
+#include <fstream>
+
+/**
+ *  This class mainly contains a function that
+ *  computes the particle relativistic kinetic energy
+ *  of each species.
+ */
+class ParticleEnergy : public ReducedDiags
+{
+public:
+
+    /** constructor
+     *  @param[in] rd_name reduced diags names */
+    ParticleEnergy(std::string rd_name);
+
+    /** This funciton computes the particle relativistic
+     *  kinetic energy (EP).
+     *  \param [in] step current time step
+     *  EP = sqrt( p^2 c^2 + m^2 c^4 ) - m c^2,
+     *  where p is the relativistic momentum,
+     *  m is the particle rest mass. */
+    virtual void ComputeDiags(int step) override final;
+
+};
+
+#endif
diff --git a/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp b/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp
new file mode 100644
index 000000000..132ad2165
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/ParticleEnergy.cpp
@@ -0,0 +1,166 @@
+/* Copyright 2019-2020 Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "ParticleEnergy.H"
+#include "WarpX.H"
+#include "WarpXConst.H"
+#include "AMReX_REAL.H"
+#include "AMReX_ParticleReduce.H"
+#include <iostream>
+#include <cmath>
+#include <limits>
+
+using namespace amrex;
+
+// constructor
+ParticleEnergy::ParticleEnergy (std::string rd_name)
+: ReducedDiags{rd_name}
+{
+    // get WarpX class object
+    auto & warpx = WarpX::GetInstance();
+
+    // get MultiParticleContainer class object
+    auto & mypc = warpx.GetPartContainer();
+
+    // get number of species (int)
+    auto nSpecies = mypc.nSpecies();
+
+    // resize data array
+    m_data.resize(2*nSpecies+2,0.0);
+
+    // get species names (std::vector<std::string>)
+    auto species_names = mypc.GetSpeciesNames();
+
+    if (ParallelDescriptor::IOProcessor())
+    {
+        if ( m_IsNotRestart )
+        {
+            // open file
+            std::ofstream ofs;
+            ofs.open(m_path + m_rd_name + "." + m_extension,
+                std::ofstream::out | std::ofstream::app);
+            // write header row
+            ofs << "#";
+            ofs << "[1]step";
+            ofs << m_sep;
+            ofs << "[2]time(s)";
+            ofs << m_sep;
+            ofs << "[3]total(J)";
+            for (int i = 0; i < nSpecies; ++i)
+            {
+                ofs << m_sep;
+                ofs << "[" + std::to_string(4+i) + "]";
+                ofs << species_names[i]+"(J)";
+            }
+            ofs << m_sep;
+            ofs << "[" + std::to_string(4+nSpecies) + "]";
+            ofs << "total.mean(J)";
+            for (int i = 0; i < nSpecies; ++i)
+            {
+                ofs << m_sep;
+                ofs << "[" + std::to_string(5+nSpecies+i) + "]";
+                ofs << species_names[i]+".mean(J)";
+            }
+            ofs << std::endl;
+            // close file
+            ofs.close();
+        }
+    }
+
+}
+// end constructor
+
+// function that computes kinetic energy
+void ParticleEnergy::ComputeDiags (int step)
+{
+
+    // Judge if the diags should be done
+    if ( (step+1) % m_freq != 0 ) { return; }
+
+    // get MultiParticleContainer class object
+    auto & mypc = WarpX::GetInstance().GetPartContainer();
+
+    // get number of species (int)
+    auto nSpecies = mypc.nSpecies();
+
+    // get species names (std::vector<std::string>)
+    auto species_names = mypc.GetSpeciesNames();
+
+    // speed of light squared
+    auto c2 = PhysConst::c * PhysConst::c;
+
+    // loop over species
+    for (int i_s = 0; i_s < nSpecies; ++i_s)
+    {
+        // get WarpXParticleContainer class object
+        auto & myspc = mypc.GetParticleContainer(i_s);
+
+        // get mass (Real)
+        auto m = myspc.getMass();
+
+        using PType = typename WarpXParticleContainer::SuperParticleType;
+
+        // Use amex::ReduceSum to compute the sum of energies of all particles
+        // held by the current MPI rank, for this species. This involves a loop over all
+        // boxes held by this MPI rank.
+        auto Etot = ReduceSum( myspc,
+        [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real
+        {
+            auto w  = p.rdata(PIdx::w);
+            auto ux = p.rdata(PIdx::ux);
+            auto uy = p.rdata(PIdx::uy);
+            auto uz = p.rdata(PIdx::uz);
+            auto us = (ux*ux + uy*uy + uz*uz);
+            return ( std::sqrt(us*c2 + c2*c2) - c2 ) * m * w;
+        });
+
+        // Same thing for the particles weights.
+        auto Wtot = ReduceSum( myspc,
+        [=] AMREX_GPU_HOST_DEVICE (const PType& p) -> Real
+        {
+            return p.rdata(PIdx::w);
+        });
+
+        // reduced sum over mpi ranks
+        ParallelDescriptor::ReduceRealSum
+            (Etot, ParallelDescriptor::IOProcessorNumber());
+        ParallelDescriptor::ReduceRealSum
+            (Wtot, ParallelDescriptor::IOProcessorNumber());
+
+        // save results for this species i_s into m_data
+        m_data[i_s+1] = Etot;
+        if ( Wtot > std::numeric_limits<Real>::min() )
+        { m_data[nSpecies+2+i_s] = Etot / Wtot; }
+        else
+        { m_data[nSpecies+2+i_s] = 0.0; }
+
+    }
+    // end loop over species
+
+    // save total energy
+    // loop over species
+    m_data[0] = 0.0;          // total energy
+    m_data[nSpecies+1] = 0.0; // total mean energy
+    for (int i_s = 0; i_s < nSpecies; ++i_s)
+    {
+        m_data[0] += m_data[i_s+1];
+        m_data[nSpecies+1] += m_data[nSpecies+2+i_s];
+    }
+    // end loop over species
+
+    /* m_data now contains up-to-date values for:
+     *  [total energy (all species),
+     *   total energy (species 1),
+     *   ...,
+     *   total energy (species n),
+     *   mean energy (all species),
+     *   mean energy (species 1),
+     *   ...,
+     *   mean energy (species n)] */
+
+}
+// end void ParticleEnergy::ComputeDiags
diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.H b/Source/Diagnostics/ReducedDiags/ReducedDiags.H
new file mode 100644
index 000000000..7ff065f49
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.H
@@ -0,0 +1,59 @@
+/* Copyright 2019-2020 Maxence Thevenet, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_REDUCEDDIAGS_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_REDUCEDDIAGS_H_
+
+#include "AMReX_REAL.H"
+#include <string>
+#include <vector>
+#include <fstream>
+
+/**
+ *  Base class for reduced diagnostics. Each type of reduced diagnostics is
+ *  implemented in a derived class, and must override the (pure virtual)
+ *  function ComputeDiags.
+ */
+class ReducedDiags
+{
+public:
+
+    /// output path (default)
+    std::string m_path = "./diags/reducedfiles/";
+
+    /// output extension (default)
+    std::string m_extension = "txt";
+
+    /// diags name
+    std::string m_rd_name;
+
+    /// output frequency
+    int m_freq = 1;
+
+    /// check if it is a restart run
+    int m_IsNotRestart = 1;
+
+    /// separator in the output file
+    std::string m_sep = ",";
+
+    /// output data
+    std::vector<amrex::Real> m_data;
+
+    /** constructor
+     *  @param[in] rd_name reduced diags name */
+    ReducedDiags(std::string rd_name);
+
+    /// function to compute diags
+    virtual void ComputeDiags(int step) = 0;
+
+    /** write to file function
+     *  @param[in] step time step */
+    virtual void WriteToFile(int step) const;
+
+};
+
+#endif
diff --git a/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp
new file mode 100644
index 000000000..81831aa79
--- /dev/null
+++ b/Source/Diagnostics/ReducedDiags/ReducedDiags.cpp
@@ -0,0 +1,92 @@
+/* Copyright 2019-2020 Maxence Thevenet, Yinjian Zhao
+ *
+ * This file is part of WarpX.
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "ReducedDiags.H"
+#include "WarpX.H"
+#include "AMReX_ParmParse.H"
+#include "AMReX_Utility.H"
+#include <iomanip>
+
+using namespace amrex;
+
+// constructor
+ReducedDiags::ReducedDiags (std::string rd_name)
+{
+
+    m_rd_name = rd_name;
+
+    ParmParse pp(m_rd_name);
+
+    // read path
+    pp.query("path", m_path);
+
+    // read extension
+    pp.query("extension", m_extension);
+
+    // creater folder
+    if (!UtilCreateDirectory(m_path, 0755))
+    { CreateDirectoryFailed(m_path); }
+
+    // check if it is a restart run
+    std::string restart_chkfile = "";
+    ParmParse pp_amr("amr");
+    pp_amr.query("restart", restart_chkfile);
+    m_IsNotRestart = restart_chkfile.empty();
+
+    // replace / create output file
+    if ( m_IsNotRestart ) // not a restart
+    {
+        std::ofstream ofs;
+        ofs.open(m_path+m_rd_name+"."+m_extension, std::ios::trunc);
+        ofs.close();
+    }
+
+    // read reduced diags frequency
+    pp.query("frequency", m_freq);
+
+    // read separator
+    pp.query("separator", m_sep);
+
+}
+// end constructor
+
+// write to file function
+void ReducedDiags::WriteToFile (int step) const
+{
+
+    // open file
+    std::ofstream ofs;
+    ofs.open(m_path + m_rd_name + "." + m_extension,
+        std::ofstream::out | std::ofstream::app);
+
+    // write step
+    ofs << step+1;
+
+    ofs << m_sep;
+
+    // set precision
+    ofs << std::fixed << std::setprecision(14) << std::scientific;
+
+    // write time
+    ofs << WarpX::GetInstance().gett_new(0);
+
+    // loop over data size and write
+    for (int i = 0; i < m_data.size(); ++i)
+    {
+        ofs << m_sep;
+        ofs << m_data[i];
+    }
+    // end loop over data size
+
+    // end line
+    ofs << std::endl;
+
+    // close file
+    ofs.close();
+
+}
+// end ReducedDiags::WriteToFile
diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp
index 84b1491c1..bb1300562 100644
--- a/Source/Evolve/WarpXEvolveEM.cpp
+++ b/Source/Evolve/WarpXEvolveEM.cpp
@@ -24,7 +24,6 @@
 #include <AMReX_AmrMeshInSituBridge.H>
 #endif
 
-
 using namespace amrex;
 
 void
@@ -206,6 +205,13 @@ WarpX::EvolveEM (int numsteps)
             t_new[i] = cur_time;
         }
 
+        /// reduced diags
+        if (reduced_diags->m_plot_rd != 0)
+        {
+            reduced_diags->ComputeDiags(step);
+            reduced_diags->WriteToFile(step);
+        }
+
         // slice gen //
         if (to_make_plot || to_write_openPMD || do_insitu || to_make_slice_plot)
         {
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index b046c667e..8b2fe1831 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -91,6 +91,13 @@ WarpX::InitData ()
 
         if ((insitu_int > 0) && (insitu_start == 0))
             UpdateInSitu();
+
+        // Write reduced diagnostics before the first iteration.
+        if (reduced_diags->m_plot_rd != 0)
+        {
+            reduced_diags->ComputeDiags(-1);
+            reduced_diags->WriteToFile(-1);
+        }
     }
 }
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index c368e2be9..ae25a8168 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -18,6 +18,7 @@
 #include <BackTransformedDiagnostic.H>
 #include <BilinearFilter.H>
 #include <NCIGodfreyFilter.H>
+#include "MultiReducedDiags.H"
 
 #ifdef WARPX_USE_PSATD
 #   include <SpectralSolver.H>
@@ -212,6 +213,9 @@ public:
     amrex::Vector<amrex::Real> mirror_z_width;
     amrex::Vector<int> mirror_z_npoints;
 
+    /// object with all reduced diagnotics, similar to MultiParticleContainer for species.
+    MultiReducedDiags* reduced_diags;
+
     void applyMirrors(amrex::Real time);
 
     void ComputeDt ();
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 71dd422ee..d93fab7df 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -201,6 +201,9 @@ WarpX::WarpX ()
     }
     do_back_transformed_particles = mypc->doBackTransformedDiagnostics();
 
+    /** create object for reduced diagnostics */
+    reduced_diags = new MultiReducedDiags();
+
     Efield_aux.resize(nlevs_max);
     Bfield_aux.resize(nlevs_max);
 
@@ -298,6 +301,8 @@ WarpX::~WarpX ()
         ClearLevel(lev);
     }
 
+    delete reduced_diags;
+
 #ifdef BL_USE_SENSEI_INSITU
     delete insitu_bridge;
 #endif