From d51a38be214cea121b608a5d2511edef0da7d4b3 Mon Sep 17 00:00:00 2001 From: Marco Garten Date: Wed, 18 Sep 2024 18:20:12 -0700 Subject: [PATCH] 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 --- .../laser_ion/CMakeLists.txt | 4 +- .../laser_ion/inputs_test_2d_laser_ion_acc | 29 ++++++--- .../inputs_test_2d_laser_ion_acc_picmi.py | 4 +- .../Physics_applications/laser_ion/plot_2d.py | 2 +- Source/Diagnostics/Diagnostics.cpp | 11 ++-- Source/Diagnostics/FullDiagnostics.H | 8 ++- Source/Diagnostics/FullDiagnostics.cpp | 65 ++++++++++++++----- 7 files changed, 86 insertions(+), 37 deletions(-) diff --git a/Examples/Physics_applications/laser_ion/CMakeLists.txt b/Examples/Physics_applications/laser_ion/CMakeLists.txt index f05203de0e8..b259421d709 100644 --- a/Examples/Physics_applications/laser_ion/CMakeLists.txt +++ b/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 ) diff --git a/Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc b/Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc index 5ad8334e9ef..6b8f77686a2 100644 --- a/Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc +++ b/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 diff --git a/Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc_picmi.py b/Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc_picmi.py index 04f9111ec5f..2cdd786aede 100755 --- a/Examples/Physics_applications/laser_ion/inputs_test_2d_laser_ion_acc_picmi.py +++ b/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, diff --git a/Examples/Physics_applications/laser_ion/plot_2d.py b/Examples/Physics_applications/laser_ion/plot_2d.py index f8a3b05d8a3..b3aefb80606 100644 --- a/Examples/Physics_applications/laser_ion/plot_2d.py +++ b/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( diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp index c8a5f93a535..feb85a6fd0b 100644 --- a/Source/Diagnostics/Diagnostics.cpp +++ b/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 ); diff --git a/Source/Diagnostics/FullDiagnostics.H b/Source/Diagnostics/FullDiagnostics.H index 9545e1314df..eebee64654a 100644 --- a/Source/Diagnostics/FullDiagnostics.H +++ b/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 diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp index 2099752a2e7..cffbc845cb0 100644 --- a/Source/Diagnostics/FullDiagnostics.cpp +++ b/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 (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,23 +260,35 @@ 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) { @@ -263,6 +296,17 @@ FullDiagnostics::DoComputeAndPack (int step, bool force_flush) 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(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 (std::round(m_average_period_time / warpx.getdt(0))); - } - - } } WARPX_ALWAYS_ASSERT_WITH_MESSAGE(