Add time step interval logic

Source/Diagnostics/Diagnostics.H

@@ -266,7 +266,11 @@ protected:
      *  The second vector is loops over the total number of levels.
      */
     amrex::Vector< amrex::Vector< amrex::MultiFab > >  m_mf_output;
-
+    /** summation multifab, where all fields (computed, cell-centered, and stacked)
+     *  are summed for every step in a time averaging period.
+     *  The first vector is for total number of snapshots. (=1 for FullDiagnostics)
+     *  The second vector is loops over the total number of levels.
+     */
     amrex::Vector< amrex::Vector <amrex::MultiFab > > m_sum_mf_output;
 
     /** Geometry that defines the domain attributes corresponding to output multifab.

Source/Diagnostics/Diagnostics.cpp

@@ -563,7 +563,8 @@ Diagnostics::ComputeAndPack ()
 
     auto & warpx = WarpX::GetInstance();
 
-    bool do_sum_now = cur_step >= time_ave 
+    // We sum up fields for later averaging and output when we are inside the averaging period
+    bool do_sum_now = cur_step >= time_ave
 
     // compute the necessary fields and store result in m_mf_output.
     for (int i_buffer = 0; i_buffer < m_num_buffers; ++i_buffer) {

Source/Diagnostics/FullDiagnostics.H

@@ -27,12 +27,16 @@ private:
      * before writing the diagnostic.
      */
     bool m_solver_deposits_current = true;
-    /** Whether the diagnostics are averaging data over time or not */
-    bool m_do_time_averaging = false;
+    /** Whether the diagnostics are averaging data over time or not
+     * Valid options are "fixed_start" and "dynamic_start".
+     */
+    std::string m_time_averaging = "none";
     /** Period to average fields over: in steps */
     unsigned int m_average_period_steps = 0;
     /** Period to average fields over: in seconds */
     amrex::Real m_average_period_seconds = 0;
+    /** Time step to start averaging */
+    unsigned int m_average_start_step = 0;
     /** Flush m_mf_output and particles to file for the i^th buffer */
     void Flush (int i_buffer, bool /* force_flush */) override;
     /** Flush raw data */

Source/Diagnostics/FullDiagnostics.cpp

@@ -142,10 +142,13 @@ FullDiagnostics::Flush ( int i_buffer, bool /* force_flush */ )
         m_output_species.at(i_buffer), nlev_output, m_file_prefix,
         m_file_min_digits, m_plot_raw_fields, m_plot_raw_fields_guards);
 
-    if (time_average) {
-       //call amrex option to 
-       m_sum_mf_output.at(i_buffer).mult(1./time_average_step);
-           m_flush_format->WriteToFile(
+    if (m_time_averaging == "fixed_start" || m_time_averaging == "dynamic_start") {
+        // Loop over the output levels and divide by the number of steps in the averaging period
+        for (int lev = 0; lev < nlev_output; ++lev) {
+            m_sum_mf_output.at(i_buffer).at(lev).mult(1./m_average_period_steps);
+        }
+
+        m_flush_format->WriteToFile(
         m_varnames, m_sum_mf_output.at(i_buffer), m_geom_output.at(i_buffer), warpx.getistep(),
         warpx.gett_new(0),
         m_output_species.at(i_buffer), nlev_output, m_file_prefix,
@@ -175,10 +178,18 @@ FullDiagnostics::DoDump (int step, int /*i_buffer*/, bool force_flush)
 bool
 FullDiagnostics::DoComputeAndPack (int step, bool force_flush)
 {
+    if (m_time_averaging == "dynamic_start")
+        m_average_start_step = m_intervals.nextContains(step) - m_average_period_steps;
+
+    // add logic here to do compute and pack if m_intervals.contains (step+1-time_average_period) or if (cur_step>time_average_startstep)
+    bool in_averaging_period = false;
+    if (step > m_intervals.nextContains(step) - m_average_start_step && step <= m_intervals.nextContains(step))
+        in_averaging_period = true;
+
     // Data must be computed and packed for full diagnostics
     // whenever the data needs to be flushed.
-    return (force_flush || m_intervals.contains(step+1));
-    // add logic here to do compute and pack if m_intervals.containes (step+1-time_average_period) or if (cur_step>time_average_startstep)
+    return (force_flush || m_intervals.contains(step+1) || in_averaging_period);
+
 }
 
 void
@@ -622,7 +633,9 @@ FullDiagnostics::InitializeBufferData (int i_buffer, int lev, bool restart ) {
     const int ngrow = (m_format == "sensei" || m_format == "ascent") ? 1 : 0;
     int const ncomp = static_cast<int>(m_varnames.size());
     m_mf_output[i_buffer][lev] = amrex::MultiFab(ba, dmap, ncomp, ngrow);
+    // Allocate MultiFab for cell-centered field output accumulation. The data will be averaged before flushing.
     m_sum_mf_output[i_buffer][lev] = amrex::MultiFab(ba, dmap, ncomp, ngrow);
+    // Initialize to zero because we add data.
     m_sum_mf_output[i_buffer][lev].setVal(0.);
     if (lev == 0) {
         // The extent of the domain covered by the diag multifab, m_mf_output

Source/Diagnostics/MultiDiagnostics.cpp

@@ -70,9 +70,10 @@ MultiDiagnostics::ReadParameters ()
         std::string diag_type_str;
         pp_diag_name.get("diag_type", diag_type_str);
         WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            diag_type_str == "Full" || diag_type_str == "BackTransformed" || diag_type_str == "BoundaryScraping",
-            "<diag>.diag_type must be Full or BackTransformed or BoundaryScraping");
+            diag_type_str == "Full" || diag_type_str == "TimeAveraged" || diag_type_str == "BackTransformed" || diag_type_str == "BoundaryScraping",
+            "<diag>.diag_type must be Full, TimeAveraged, BackTransformed or BoundaryScraping");
         if (diag_type_str == "Full") { diags_types[i] = DiagTypes::Full; }
+        if (diag_type_str == "TimeAveraged") { diags_types[i] = DiagTypes::TimeAveraged; }
         if (diag_type_str == "BackTransformed") { diags_types[i] = DiagTypes::BackTransformed; }
         if (diag_type_str == "BoundaryScraping") { diags_types[i] = DiagTypes::BoundaryScraping; }
     }