diff --git a/Source/Diagnostics/Diagnostics.H b/Source/Diagnostics/Diagnostics.H index 81553433dde..68575bbd42c 100644 --- a/Source/Diagnostics/Diagnostics.H +++ b/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 > m_sum_mf_output; /** Geometry that defines the domain attributes corresponding to output multifab. diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index e4e4a98f032..ae40044adf9 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/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) { diff --git a/Source/Diagnostics/FullDiagnostics.H b/Source/Diagnostics/FullDiagnostics.H index a58c2485c8f..e5788dbf20a 100644 --- a/Source/Diagnostics/FullDiagnostics.H +++ b/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 */ diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index 9e207a03bd9..da6b8523491 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/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(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 diff --git a/Source/Diagnostics/MultiDiagnostics.cpp b/Source/Diagnostics/MultiDiagnostics.cpp index ecba8831b20..739d5ed330c 100644 --- a/Source/Diagnostics/MultiDiagnostics.cpp +++ b/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_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_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; } }