Add time averaged diags to laser ion example

- added time averaged diags to laser-ion acceleration example
- fix first issues that came up when testing this

Examples/Physics_applications/laser_ion/CMakeLists.txt

@@ -7,7 +7,7 @@ add_warpx_test(
     2  # nprocs
     inputs_test_2d_laser_ion_acc  # inputs
     analysis_default_openpmd_regression.py  # analysis
-    diags/diag1/  # output
+    diags/diagInst/  # output
     OFF  # dependency
 )
 
@@ -17,6 +17,6 @@ add_warpx_test(
     2  # nprocs
     inputs_test_2d_laser_ion_acc_picmi.py  # inputs
     analysis_default_openpmd_regression.py  # analysis
-    diags/diag1/  # output
+    diags/diagInst/  # output
     OFF  # dependency
 )

Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc

@@ -200,18 +200,29 @@ laser1.profile_focal_distance = 4.0e-6  # focal distance from the antenna [m]
 #################################
 # Diagnostics
 #
-diagnostics.diags_names = diag1 openPMDfw openPMDbw
+diagnostics.diags_names = diagInst diagTimeAvg openPMDfw openPMDbw
 
-diag1.intervals = 100
-diag1.diag_type = Full
-diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen
+diagInst.intervals = 100
+diagInst.diag_type = Full
+diagInst.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen
 # reduce resolution of output fields
-diag1.coarsening_ratio = 4 4
+diagInst.coarsening_ratio = 4 4
 # demonstration of a spatial and momentum filter
-diag1.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)
-diag1.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)
-diag1.format = openpmd
-diag1.openpmd_backend = h5
+diagInst.electrons.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)
+diagInst.hydrogen.plot_filter_function(t,x,y,z,ux,uy,uz) = (uz>=0) * (x<1.0e-6) * (x>-1.0e-6)
+diagInst.format = openpmd
+diagInst.openpmd_backend = h5
+
+diagTimeAvg.intervals = 100::100  # TODO just write 100 and make step 0 an instantaneous diagnostic and do not forget to document
+diagTimeAvg.diag_type = TimeAveraged
+diagTimeAvg.time_average_mode = dynamic_start
+diagTimeAvg.average_period_time = 2.67e-15  # period of 800 nm light waves
+diagTimeAvg.write_species = 0
+diagTimeAvg.fields_to_plot = Ey Ez Bx By Bz jx jy jz rho rho_electrons rho_hydrogen
+# reduce resolution of output fields
+diagTimeAvg.coarsening_ratio = 4 4
+diagTimeAvg.format = openpmd
+diagTimeAvg.openpmd_backend = h5
 
 openPMDfw.intervals = 100
 openPMDfw.diag_type = Full

Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc_picmi.py

@@ -140,7 +140,7 @@
 
 # Diagnostics
 particle_diag = picmi.ParticleDiagnostic(
-    name="diag1",
+    name="diagInst",
     period=100,
     warpx_format="openpmd",
     warpx_openpmd_backend="h5",
@@ -153,7 +153,7 @@
 for ncell_comp, cr in zip([nx, nz], coarsening_ratio):
     ncell_field.append(int(ncell_comp / cr))
 field_diag = picmi.FieldDiagnostic(
-    name="diag1",
+    name="diagInst",
     grid=grid,
     period=100,
     number_of_cells=ncell_field,

Examples/Physics_applications/laser_ion/plot_2d.py

@@ -259,7 +259,7 @@ def visualize_particle_histogram_iteration(
         "-d",
         "--diag_dir",
         type=str,
-        default="./diags/diag1",
+        default="./diags/diagInst",
         help="Directory containing density and field diagnostics",
     )
     parser.add_argument(

Source/Diagnostics/Diagnostics.cpp

@@ -534,11 +534,14 @@ Diagnostics::InitBaseData ()
     for (int i = 0; i < m_num_buffers; ++i) {
         m_mf_output[i].resize( nmax_lev );
     }
-    m_sum_mf_output.resize(m_num_buffers);
-    for (int i = 0; i < m_num_buffers; ++i) {
-        m_sum_mf_output[i].resize( nmax_lev );
-    }
 
+    // allocate vector of buffers and vector of levels for each buffer for summation multifab for TimeAveragedDiagnostics
+    if (m_diag_type == DiagTypes::TimeAveraged) {
+        m_sum_mf_output.resize(m_num_buffers);
+        for (int i = 0; i < m_num_buffers; ++i) {
+            m_sum_mf_output[i].resize( nmax_lev );
+        }
+    }
 
     // allocate vector of geometry objects corresponding to each output multifab.
     m_geom_output.resize( m_num_buffers );

Source/Diagnostics/FullDiagnostics.H

@@ -38,15 +38,17 @@ private:
     /** Whether the diagnostics are averaging data over time or not */
     TimeAverageType m_time_average_type = TimeAverageType::None;
     /** Period to average fields over: in steps */
-    int m_average_period_steps = 0;
+    int m_average_period_steps = -1;
     /** Period to average fields over: in seconds */
-    amrex::Real m_average_period_time = 0;
+    amrex::Real m_average_period_time = -1.0;
     /** Time step to start averaging */
-    int m_average_start_step = 0;
+    int m_average_start_step = -1;
     /** Flush m_mf_output or m_sum_mf_output and particles to file for the i^th buffer */
     void Flush (int i_buffer, bool /* force_flush */) override;
     /** Flush raw data */
     void FlushRaw ();
+    /** Initialize Data required to compute TimeAveraged diagnostics */
+    void DerivedInitData () override;
     /** whether to compute and pack cell-centered data in m_mf_output
      * \param[in] step current time step
      * \param[in] force_flush if true, return true for any step since output must be

Source/Diagnostics/FullDiagnostics.cpp

@@ -55,6 +55,27 @@ FullDiagnostics::FullDiagnostics (int i, const std::string& name, DiagTypes diag
     BackwardCompatibility();
 }
 
+void
+FullDiagnostics::DerivedInitData() {
+    if (m_diag_type == DiagTypes::TimeAveraged) {
+        auto & warpx = WarpX::GetInstance();
+        if (m_time_average_type == TimeAverageType::Dynamic) {
+
+            // already checked in ReadParameters that only one of the parameters is set
+            // calculate the other averaging period parameter from the other one, respectively
+            if (m_average_period_steps > 0) {
+                m_average_period_time = m_average_period_steps * warpx.getdt(0);
+            } else if (m_average_period_time > 0) {
+                m_average_period_steps = static_cast<int> (std::round(m_average_period_time / warpx.getdt(0)));
+            }
+            amrex::Print() << Utils::TextMsg::Info(
+                "Initializing TimeAveragedDiagnostics " + m_diag_name
+                + " with an averaging period of " + std::to_string(m_average_period_steps) + " steps"
+                );
+        }
+    }
+}
+
 void
 FullDiagnostics::InitializeParticleBuffer ()
 {
@@ -239,30 +260,53 @@ FullDiagnostics::DoComputeAndPack (int step, bool force_flush)
     // Start averaging at output step (from diag.intervals) - period + 1
     bool in_averaging_period = false;
     if (m_diag_type == DiagTypes::TimeAveraged) {
+
+        // TODO we need a good solution to avoid that there is output at step 0
+        if (step == 0) {
+            return false;
+        }
+
         if (m_time_average_type == TimeAverageType::Dynamic) {
             m_average_start_step = m_intervals.nextContains(step) - m_average_period_steps;
             // check that the periods do not overlap and that the start step is not negative
-            if (m_average_start_step > 0) {
-                if (m_average_start_step < m_intervals.previousContains(step)) {
+            if (m_average_start_step >= 0) {
+                // The start step cannot be on an interval step because then we would begin a new period and also output the old one
+                if (m_average_start_step <= m_intervals.previousContains(step)) {
                     WARPX_ABORT_WITH_MESSAGE(
                             "Averaging periods may not overlap within a single diagnostic. "
                             "Please create a second diagnostic for overlapping time averaging periods "
                             "and account for the increased memory consumption."
                     );
-                } else {
-                    WARPX_ABORT_WITH_MESSAGE("The step to begin time averaging may not be a negative number.");
                 }
+
+            } else {
+                WARPX_ABORT_WITH_MESSAGE(
+                        "The step to begin time averaging ("
+                        + std::to_string(m_average_start_step)
+                        + ") may not be a negative number."
+                        );
             }
         }
 
-        if (step > m_intervals.nextContains(step) - m_average_start_step && step <= m_intervals.nextContains(step)) {
+        if (step >= m_average_start_step && step <= m_intervals.nextContains(step)) {
             in_averaging_period = true;
 
             if (m_time_average_type == TimeAverageType::Static) {
                 // Update time averaging period to current step
                 m_average_period_steps = step - m_average_start_step;
             }
         }
+        // Print information on when time-averaging is active
+        if (in_averaging_period) {
+            if (step == m_average_start_step) {
+                amrex::Print() << Utils::TextMsg::Info(
+                        "Begin time averaging for " + m_diag_name + " and output at step "
+                        + std::to_string(m_intervals.nextContains(step))
+                );
+            } else {
+                amrex::Print() << Utils::TextMsg::Info("Time-averaging during this step for diagnostic: " + m_diag_name);
+            }
+        }
     }
 
     // Data must be computed and packed for full diagnostics
@@ -700,17 +744,6 @@ FullDiagnostics::InitializeBufferData (int i_buffer, int lev, bool restart ) {
                 (ba.getCellCenteredBox( static_cast<int>(ba.size())-1 ).bigEnd(idim) + 1) * warpx.Geom(lev).CellSize(idim));
         }
 
-        if (m_time_average_type == TimeAverageType::Dynamic) {
-
-            // already checked in ReadParameters that only one of the parameters is set
-            // calculate the other averaging period parameter from the other one, respectively
-            if (m_average_period_steps > 0) {
-                m_average_period_time = m_average_period_steps * warpx.getdt(0);
-            } else if (m_average_period_time > 0) {
-                m_average_period_steps = static_cast<int> (std::round(m_average_period_time / warpx.getdt(0)));
-            }
-
-        }
     }
 
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(