From 22ed23940fb66426e71e29b5ab80fd3b45a12401 Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Wed, 8 Jan 2025 19:32:49 -0800 Subject: [PATCH 1/8] add diff lumi diag 2d --- .../Diagnostics/ReducedDiags/CMakeLists.txt | 1 + .../ReducedDiags/DifferentialLuminosity2D.H | 65 +++ .../ReducedDiags/DifferentialLuminosity2D.cpp | 416 ++++++++++++++++++ Source/Diagnostics/ReducedDiags/Make.package | 1 + .../ReducedDiags/MultiReducedDiags.cpp | 40 +- 5 files changed, 504 insertions(+), 19 deletions(-) create mode 100644 Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H create mode 100644 Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp diff --git a/Source/Diagnostics/ReducedDiags/CMakeLists.txt b/Source/Diagnostics/ReducedDiags/CMakeLists.txt index bbf1b6b65b0..9514c860945 100644 --- a/Source/Diagnostics/ReducedDiags/CMakeLists.txt +++ b/Source/Diagnostics/ReducedDiags/CMakeLists.txt @@ -6,6 +6,7 @@ foreach(D IN LISTS WarpX_DIMS) ChargeOnEB.cpp ColliderRelevant.cpp DifferentialLuminosity.cpp + DifferentialLuminosity2D.cpp FieldEnergy.cpp FieldMaximum.cpp FieldMomentum.cpp diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H new file mode 100644 index 00000000000..865169ee380 --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H @@ -0,0 +1,65 @@ +/* Copyright 2023 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Arianna Formenti, Remi Lehe + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_ +#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_ + +#include "ReducedDiags.H" +#include +#include + +#include +#include +#include + +/** + * This class contains the differential luminosity diagnostics. + */ +class DifferentialLuminosity2D : public ReducedDiags +{ +public: + + /** + * constructor + * @param[in] rd_name reduced diags names + */ + DifferentialLuminosity2D(const std::string& rd_name); + + /// File type + std::string m_openpmd_backend {"default"}; + + /// minimum number of digits for file suffix (file-based only supported for now) */ + int m_file_min_digits = 6; + + /// name of the two colliding species + std::vector m_beam_name; + + /// number of bins for the c.o.m. energy of the 2 species + int m_bin_num_1; + int m_bin_num_2; + + /// max and min bin values + amrex::Real m_bin_max_1; + amrex::Real m_bin_min_1; + amrex::Real m_bin_max_2; + amrex::Real m_bin_min_2; + + /// bin size + amrex::Real m_bin_size_1; + amrex::Real m_bin_size_2; + + /// output data + amrex::TableData m_h_data_2D; + + void ComputeDiags(int step) final; + + void WriteToFile (int step) const final; + +}; + +#endif // WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_ diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp new file mode 100644 index 00000000000..e90bcaac77c --- /dev/null +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp @@ -0,0 +1,416 @@ +/* Copyright 2023 The WarpX Community + * + * This file is part of WarpX. + * + * Authors: Arianna Formenti, Yinjian Zhao, Remi Lehe + * License: BSD-3-Clause-LBNL + */ +#include "DifferentialLuminosity2D.H" + +#include "Diagnostics/ReducedDiags/ReducedDiags.H" +#include "Diagnostics/OpenPMDHelpFunction.H" +#include "Particles/MultiParticleContainer.H" +#include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/SpeciesPhysicalProperties.H" +#include "Particles/WarpXParticleContainer.H" +#include "Utils/ParticleUtils.H" +#include "Utils/Parser/ParserUtils.H" +#include "Utils/WarpXConst.H" +#include "Utils/TextMsg.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef WARPX_USE_OPENPMD +# include +#endif + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using ParticleType = WarpXParticleContainer::ParticleType; +using ParticleTileType = WarpXParticleContainer::ParticleTileType; +using ParticleTileDataType = ParticleTileType::ParticleTileDataType; +using ParticleBins = amrex::DenseBins; +using index_type = ParticleBins::index_type; + +#ifdef WARPX_USE_OPENPMD +namespace io = openPMD; +#endif + +using namespace amrex; + +DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) +: ReducedDiags{rd_name} +{ + // RZ coordinate is not supported +#if (defined WARPX_DIM_RZ) + WARPX_ABORT_WITH_MESSAGE( + "DifferentialLuminosity2D diagnostics do not work in RZ geometry."); +#endif + + // read colliding species names - must be 2 + amrex::ParmParse pp_rd_name(m_rd_name); + pp_rd_name.getarr("species", m_beam_name); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + m_beam_name.size() == 2u, + "DifferentialLuminosity2D diagnostics must involve exactly two species"); + + // should we add an assert that checks that the selected species actually exist? + // ref. ParticleHistrogram2D + + pp_rd_name.query("openpmd_backend", m_openpmd_backend); + pp_rd_name.query("file_min_digits", m_file_min_digits); + // pick first available backend if default is chosen + if( m_openpmd_backend == "default" ) { + m_openpmd_backend = WarpXOpenPMDFileType(); + } + pp_rd_name.add("openpmd_backend", m_openpmd_backend); + + // read bin parameters for first species + int bin_num_1 = 0; + amrex::Real bin_max_1 = 0.0_rt, bin_min_1 = 0.0_rt; + utils::parser::getWithParser(pp_rd_name, "bin_number_1", bin_num_1); + utils::parser::getWithParser(pp_rd_name, "bin_max_1", bin_max_1); + utils::parser::getWithParser(pp_rd_name, "bin_min_1", bin_min_1); + m_bin_num_1 = bin_num_1; + m_bin_max_1 = bin_max_1; + m_bin_min_1 = bin_min_1; + m_bin_size_1 = (bin_max_1 - bin_min_1) / bin_num_1; + + // read bin parameters for second species + int bin_num_2 = 0; + amrex::Real bin_max_2 = 0.0_rt, bin_min_2 = 0.0_rt; + utils::parser::getWithParser(pp_rd_name, "bin_number_2", bin_num_2); + utils::parser::getWithParser(pp_rd_name, "bin_max_2", bin_max_2); + utils::parser::getWithParser(pp_rd_name, "bin_min_2", bin_min_2); + m_bin_num_2 = bin_num_2; + m_bin_max_2 = bin_max_2; + m_bin_min_2 = bin_min_2; + m_bin_size_2 = (bin_max_2 - bin_min_2) / bin_num_2; + +} + +void DifferentialLuminosity2D::ComputeDiags (int step) +{ +#if defined(WARPX_DIM_RZ) + amrex::ignore_unused(step); +#else + + // resize data array + Array tlo{0,0}; // lower bounds + Array thi{m_bin_num_1-1, m_bin_num_2-1}; // inclusive upper bounds + amrex::TableData d_data_2D(tlo, thi); + m_h_data_2D.resize(tlo, thi, The_Pinned_Arena()); + auto const& h_table_data = m_h_data_2D.table(); + + // Initialize data on the host + if(step==0){ + for (int i = tlo[0]; i <= thi[0]; ++i) { + for (int j = tlo[1]; j <= thi[1]; ++j) { + h_table_data(i,j) = 0.0_rt; + } + }} + + d_data_2D.copy(m_h_data_2D); + auto d_table = d_data_2D.table(); + + WARPX_PROFILE("DifferentialLuminosity2D::ComputeDiags"); + + // Since this diagnostic *accumulates* the luminosity in the + // array d_data, we add contributions at *each timestep*, but + // we only write the data to file at intervals specified by the user. + const Real c_sq = PhysConst::c*PhysConst::c; + const Real c_over_qe = PhysConst::c/PhysConst::q_e; + + // get a reference to WarpX instance + auto& warpx = WarpX::GetInstance(); + const Real dt = warpx.getdt(0); + // get cell volume + Geometry const & geom = warpx.Geom(0); + const Real dV = AMREX_D_TERM(geom.CellSize(0), *geom.CellSize(1), *geom.CellSize(2)); + + // declare local variables + auto const num_bins_1 = m_bin_num_1; + Real const bin_min_1 = m_bin_min_1; + Real const bin_size_1 = m_bin_size_1; + auto const num_bins_2 = m_bin_num_2; + Real const bin_min_2 = m_bin_min_2; + Real const bin_size_2 = m_bin_size_2; + + // get MultiParticleContainer class object + const MultiParticleContainer& mypc = warpx.GetPartContainer(); + + auto& species_1 = mypc.GetParticleContainerFromName(m_beam_name[0]); + auto& species_2 = mypc.GetParticleContainerFromName(m_beam_name[1]); + + const ParticleReal m1 = species_1.getMass(); + const ParticleReal m2 = species_2.getMass(); + + // Enable tiling + amrex::MFItInfo info; + if (amrex::Gpu::notInLaunchRegion()) { info.EnableTiling(WarpXParticleContainer::tile_size); } + + int const nlevs = std::max(0, species_1.finestLevel()+1); // species_1 ? + for (int lev = 0; lev < nlevs; ++lev) { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + + for (amrex::MFIter mfi = species_1.MakeMFIter(lev, info); mfi.isValid(); ++mfi){ + + ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi); + ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi); + + ParticleBins bins_1 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_1 ); + ParticleBins bins_2 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_2 ); + + // Species + const auto soa_1 = ptile_1.getParticleTileData(); + index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr(); + index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr(); + + // Particle data in the tile/box + amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; + amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; + amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; // v*gamma=p/m + amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; + bool const species1_is_photon = species_1.AmIA(); + + const auto soa_2 = ptile_2.getParticleTileData(); + index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr(); + index_type const* AMREX_RESTRICT cell_offsets_2 = bins_2.offsetsPtr(); + + amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; + amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux]; + amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy]; + amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; + bool const species2_is_photon = species_2.AmIA(); + + // Extract low-level data + auto const n_cells = static_cast(bins_1.numBins()); + + // Loop over cells + amrex::ParallelFor( n_cells, + [=] AMREX_GPU_DEVICE (int i_cell) noexcept + { + + // The particles from species1 that are in the cell `i_cell` are + // given by the `indices_1[cell_start_1:cell_stop_1]` + index_type const cell_start_1 = cell_offsets_1[i_cell]; + index_type const cell_stop_1 = cell_offsets_1[i_cell+1]; + // Same for species 2 + index_type const cell_start_2 = cell_offsets_2[i_cell]; + index_type const cell_stop_2 = cell_offsets_2[i_cell+1]; + + for(index_type i_1=cell_start_1; i_1=num_bins_1 ) { continue; } // discard if out-of-range + + // determine particle 2 bin + int const bin_2 = int(Math::floor((E_2-bin_min_2)/bin_size_2)); + if ( bin_2<0 || bin_2>=num_bins_2 ) { continue; } // discard if out-of-range + + + Real const inv_p1t = 1.0_rt/p1t; + Real const inv_p2t = 1.0_rt/p2t; + + Real const beta1_sq = (p1x*p1x + p1y*p1y + p1z*p1z) * inv_p1t*inv_p1t; + Real const beta2_sq = (p2x*p2x + p2y*p2y + p2z*p2z) * inv_p2t*inv_p2t; + Real const beta1_dot_beta2 = (p1x*p2x + p1y*p2y + p1z*p2z) * inv_p1t*inv_p2t; + + // Here we use the fact that: + // (v1 - v2)^2 = v1^2 + v2^2 - 2 v1.v2 + // and (v1 x v2)^2 = v1^2 v2^2 - (v1.v2)^2 + // we also use beta=v/c instead of v + + Real const radicand = beta1_sq + beta2_sq - 2*beta1_dot_beta2 - beta1_sq*beta2_sq + beta1_dot_beta2*beta1_dot_beta2; + + Real const d2L_dE1_dE2 = PhysConst::c * std::sqrt( radicand ) * w1[j_1] * w2[j_2] / (dV * bin_size_1 * bin_size_2) * dt; // m^-2 eV^-2 + if (d2L_dE1_dE2 < 0){ + amrex::Print(0) << "NEGATIVE VALUE!!!!!!!!!!!!!!!!!!!!!!!!!\n"; + } + amrex::Real &data = d_table(bin_1, bin_2); + amrex::HostDevice::Atomic::Add(&data, d2L_dE1_dE2); + + } // particles species 2 + } // particles species 1 + }); // cells + } // boxes + } // levels + + + // Only write to file at intervals specified by the user. + // At these intervals, the data needs to ready on the CPU, + // so we copy it from the GPU to the CPU and reduce across MPI ranks. + if (m_intervals.contains(step+1)) { + + // Copy data from GPU memory + m_h_data_2D.copy(d_data_2D); + + // reduced sum over mpi ranks + const int size = static_cast (d_data_2D.size()); + ParallelDescriptor::ReduceRealSum + (h_table_data.p, size, ParallelDescriptor::IOProcessorNumber()); + } + + // Return for all that are not IO processor + if ( !ParallelDescriptor::IOProcessor() ) { return; } + +#endif // not RZ +} // end void DifferentialLuminosity2D::ComputeDiags + + + + +void DifferentialLuminosity2D::WriteToFile (int step) const +{ + // Judge if the diags should be done at this step + if (!m_intervals.contains(step+1)) { return; } + +#ifdef WARPX_USE_OPENPMD + // only IO processor writes + if ( !ParallelDescriptor::IOProcessor() ) { return; } + + // TODO: support different filename templates + std::string filename = "openpmd"; + // TODO: support also group-based encoding + const std::string fileSuffix = std::string("_%0") + std::to_string(m_file_min_digits) + std::string("T"); + filename = filename.append(fileSuffix).append(".").append(m_openpmd_backend); + + // transform paths for Windows + #ifdef _WIN32 + const std::string filepath = openPMD::auxiliary::replace_all( + m_path + m_rd_name + "/" + filename, "/", "\\"); + #else + const std::string filepath = m_path + m_rd_name + "/" + filename; + #endif + + // Create the OpenPMD series + auto series = io::Series( + filepath, + io::Access::CREATE); + auto i = series.iterations[step + 1]; + // record + auto f_mesh = i.meshes["data"]; + //f_mesh.setAttribute("energy_1", "E1"); + //f_mesh.setAttribute("energy_2", "E2"); + f_mesh.setUnitDimension({ + {io::UnitDimension::L, -6}, + {io::UnitDimension::M, -2}, + {io::UnitDimension::T, 4} + }); + + // add units! + + // record components + auto data = f_mesh[io::RecordComponent::SCALAR]; + + data.setUnitSI(std::pow(PhysConst::q_e, -2)); + + // meta data + f_mesh.setAxisLabels({"E2", "E1"}); // ordinate, abscissa + std::vector< double > const& gridGlobalOffset = {m_bin_min_2, m_bin_min_1}; + f_mesh.setGridGlobalOffset(gridGlobalOffset); + f_mesh.setGridSpacing({m_bin_size_2, m_bin_size_1}); + + data.setPosition({0.5, 0.5}); + + auto dataset = io::Dataset( + io::determineDatatype(), + {static_cast(m_bin_num_2), static_cast(m_bin_num_1)}); + data.resetDataset(dataset); + + // UNIT DIMENSION IS NOT SET ON THE VALUES + + // Get time at level 0 + auto & warpx = WarpX::GetInstance(); + auto const time = warpx.gett_new(0); + i.setTime(time); + + auto const& h_table_data = m_h_data_2D.table(); + data.storeChunkRaw( + h_table_data.p, + {0, 0}, + {static_cast(m_bin_num_2), static_cast(m_bin_num_1)}); + + series.flush(); + i.close(); + series.close(); +#else + amrex::ignore_unused(step); + WARPX_ABORT_WITH_MESSAGE("ParticleHistogram2D: Needs openPMD-api compiled into WarpX, but was not found!"); +#endif +} diff --git a/Source/Diagnostics/ReducedDiags/Make.package b/Source/Diagnostics/ReducedDiags/Make.package index 2611831a3dd..335b2396d5b 100644 --- a/Source/Diagnostics/ReducedDiags/Make.package +++ b/Source/Diagnostics/ReducedDiags/Make.package @@ -4,6 +4,7 @@ CEXE_sources += BeamRelevant.cpp CEXE_sources += ChargeOnEB.cpp CEXE_sources += ColliderRelevant.cpp CEXE_sources += DifferentialLuminosity.cpp +CEXE_sources += DifferentialLuminosity2D.cpp CEXE_sources += FieldEnergy.cpp CEXE_sources += FieldMaximum.cpp CEXE_sources += FieldMomentum.cpp diff --git a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp index 5035eac58a8..2239aef8028 100644 --- a/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp +++ b/Source/Diagnostics/ReducedDiags/MultiReducedDiags.cpp @@ -10,6 +10,7 @@ #include "ChargeOnEB.H" #include "ColliderRelevant.H" #include "DifferentialLuminosity.H" +#include "DifferentialLuminosity2D.H" #include "FieldEnergy.H" #include "FieldMaximum.H" #include "FieldMomentum.H" @@ -53,25 +54,26 @@ MultiReducedDiags::MultiReducedDiags () using CS = const std::string& ; const auto reduced_diags_dictionary = std::map(CS)>>{ - {"BeamRelevant", [](CS s){return std::make_unique(s);}}, - {"ChargeOnEB", [](CS s){return std::make_unique(s);}}, - {"ColliderRelevant", [](CS s){return std::make_unique(s);}}, - {"DifferentialLuminosity",[](CS s){return std::make_unique(s);}}, - {"ParticleEnergy", [](CS s){return std::make_unique(s);}}, - {"ParticleExtrema", [](CS s){return std::make_unique(s);}}, - {"ParticleHistogram", [](CS s){return std::make_unique(s);}}, - {"ParticleHistogram2D", [](CS s){return std::make_unique(s);}}, - {"ParticleMomentum", [](CS s){return std::make_unique(s);}}, - {"ParticleNumber", [](CS s){return std::make_unique(s);}}, - {"FieldEnergy", [](CS s){return std::make_unique(s);}}, - {"FieldMaximum", [](CS s){return std::make_unique(s);}}, - {"FieldMomentum", [](CS s){return std::make_unique(s);}}, - {"FieldProbe", [](CS s){return std::make_unique(s);}}, - {"FieldReduction", [](CS s){return std::make_unique(s);}}, - {"LoadBalanceCosts", [](CS s){return std::make_unique(s);}}, - {"LoadBalanceEfficiency", [](CS s){return std::make_unique(s);}}, - {"RhoMaximum", [](CS s){return std::make_unique(s);}}, - {"Timestep", [](CS s){return std::make_unique(s);}} + {"BeamRelevant", [](CS s){return std::make_unique(s);}}, + {"ChargeOnEB", [](CS s){return std::make_unique(s);}}, + {"ColliderRelevant", [](CS s){return std::make_unique(s);}}, + {"DifferentialLuminosity", [](CS s){return std::make_unique(s);}}, + {"DifferentialLuminosity2D",[](CS s){return std::make_unique(s);}}, + {"ParticleEnergy", [](CS s){return std::make_unique(s);}}, + {"ParticleExtrema", [](CS s){return std::make_unique(s);}}, + {"ParticleHistogram", [](CS s){return std::make_unique(s);}}, + {"ParticleHistogram2D", [](CS s){return std::make_unique(s);}}, + {"ParticleMomentum", [](CS s){return std::make_unique(s);}}, + {"ParticleNumber", [](CS s){return std::make_unique(s);}}, + {"FieldEnergy", [](CS s){return std::make_unique(s);}}, + {"FieldMaximum", [](CS s){return std::make_unique(s);}}, + {"FieldMomentum", [](CS s){return std::make_unique(s);}}, + {"FieldProbe", [](CS s){return std::make_unique(s);}}, + {"FieldReduction", [](CS s){return std::make_unique(s);}}, + {"LoadBalanceCosts", [](CS s){return std::make_unique(s);}}, + {"LoadBalanceEfficiency", [](CS s){return std::make_unique(s);}}, + {"RhoMaximum", [](CS s){return std::make_unique(s);}}, + {"Timestep", [](CS s){return std::make_unique(s);}} }; // loop over all reduced diags and fill m_multi_rd with requested reduced diags std::transform(m_rd_names.begin(), m_rd_names.end(), std::back_inserter(m_multi_rd), From d6b92de93bf0bc676d6b47fae280e900bb57554f Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Wed, 8 Jan 2025 19:33:17 -0800 Subject: [PATCH 2/8] add test --- Examples/Tests/diff_lumi_diag/CMakeLists.txt | 5 ++++- Examples/Tests/diff_lumi_diag/inputs_base_3d | 12 +++++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/Examples/Tests/diff_lumi_diag/CMakeLists.txt b/Examples/Tests/diff_lumi_diag/CMakeLists.txt index f16449a976c..9a4e58d0e62 100644 --- a/Examples/Tests/diff_lumi_diag/CMakeLists.txt +++ b/Examples/Tests/diff_lumi_diag/CMakeLists.txt @@ -1,6 +1,6 @@ # Add tests (alphabetical order) ############################################## # - +if(WarpX_FFT) add_warpx_test( test_3d_diff_lumi_diag_leptons # name 3 # dims @@ -10,7 +10,9 @@ add_warpx_test( "analysis_default_regression.py --path diags/diag1000080 --rtol 1e-2" # checksum OFF # dependency ) +endif() +if(WarpX_FFT) add_warpx_test( test_3d_diff_lumi_diag_photons # name 3 # dims @@ -20,3 +22,4 @@ add_warpx_test( "analysis_default_regression.py --path diags/diag1000080 --rtol 1e-2" # checksum OFF # dependency ) +endif() diff --git a/Examples/Tests/diff_lumi_diag/inputs_base_3d b/Examples/Tests/diff_lumi_diag/inputs_base_3d index ba3c823b52b..4724fec3af0 100644 --- a/Examples/Tests/diff_lumi_diag/inputs_base_3d +++ b/Examples/Tests/diff_lumi_diag/inputs_base_3d @@ -93,7 +93,7 @@ diag1.dump_last_timestep = 1 diag1.species = beam1 beam2 # REDUCED -warpx.reduced_diags_names = DifferentialLuminosity_beam1_beam2 +warpx.reduced_diags_names = DifferentialLuminosity_beam1_beam2 DifferentialLuminosity2d_beam1_beam2 DifferentialLuminosity_beam1_beam2.type = DifferentialLuminosity DifferentialLuminosity_beam1_beam2.intervals = 5 @@ -101,3 +101,13 @@ DifferentialLuminosity_beam1_beam2.species = beam1 beam2 DifferentialLuminosity_beam1_beam2.bin_number = 128 DifferentialLuminosity_beam1_beam2.bin_max = 2.1*beam_energy_eV DifferentialLuminosity_beam1_beam2.bin_min = 1.9*beam_energy_eV + +DifferentialLuminosity2d_beam1_beam2.type = DifferentialLuminosity2D +DifferentialLuminosity2d_beam1_beam2.intervals = 5 +DifferentialLuminosity2d_beam1_beam2.species = beam1 beam2 +DifferentialLuminosity2d_beam1_beam2.bin_number_1 = 128 +DifferentialLuminosity2d_beam1_beam2.bin_max_1 = 1.1*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_min_1 = 0.9*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_number_2 = 256 +DifferentialLuminosity2d_beam1_beam2.bin_max_2 = 1.1*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_min_2 = 0.9*beam_energy_eV From 011597164941e98e0333ea8e30ac2c307003e6dc Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Thu, 9 Jan 2025 18:49:21 -0800 Subject: [PATCH 3/8] moved device vector, rm openpmd setUnitsSI --- .../ReducedDiags/DifferentialLuminosity2D.H | 5 +++ .../ReducedDiags/DifferentialLuminosity2D.cpp | 44 ++++++++++--------- 2 files changed, 29 insertions(+), 20 deletions(-) diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H index 865169ee380..d3c776ee9f8 100644 --- a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H @@ -60,6 +60,11 @@ public: void WriteToFile (int step) const final; +private: + + /// output data // Array in which to accumulate the luminosity across timesteps + amrex::TableData m_d_data_2D; + }; #endif // WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_ diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp index e90bcaac77c..b66e036934c 100644 --- a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp @@ -122,34 +122,38 @@ DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) m_bin_min_2 = bin_min_2; m_bin_size_2 = (bin_max_2 - bin_min_2) / bin_num_2; -} - -void DifferentialLuminosity2D::ComputeDiags (int step) -{ -#if defined(WARPX_DIM_RZ) - amrex::ignore_unused(step); -#else // resize data array Array tlo{0,0}; // lower bounds Array thi{m_bin_num_1-1, m_bin_num_2-1}; // inclusive upper bounds - amrex::TableData d_data_2D(tlo, thi); m_h_data_2D.resize(tlo, thi, The_Pinned_Arena()); + + auto const& h_table_data = m_h_data_2D.table(); // Initialize data on the host - if(step==0){ for (int i = tlo[0]; i <= thi[0]; ++i) { for (int j = tlo[1]; j <= thi[1]; ++j) { h_table_data(i,j) = 0.0_rt; } - }} + } + + m_d_data_2D.resize(tlo, thi); + m_d_data_2D.copy(m_h_data_2D); + Gpu::streamSynchronize(); - d_data_2D.copy(m_h_data_2D); - auto d_table = d_data_2D.table(); +} + +void DifferentialLuminosity2D::ComputeDiags (int step) +{ +#if defined(WARPX_DIM_RZ) + amrex::ignore_unused(step); +#else WARPX_PROFILE("DifferentialLuminosity2D::ComputeDiags"); + auto d_table = m_d_data_2D.table(); + // Since this diagnostic *accumulates* the luminosity in the // array d_data, we add contributions at *each timestep*, but // we only write the data to file at intervals specified by the user. @@ -261,7 +265,7 @@ void DifferentialLuminosity2D::ComputeDiags (int step) Real const u2_sq = u2x[j_2]*u2x[j_2] + u2y[j_2]*u2y[j_2] + u2z[j_2]*u2z[j_2]; if (species2_is_photon) { // photon case (momentum is normalized by m_e in WarpX) - p2t = PhysConst::m_e*std::sqrt(u2_sq); + p2t = PhysConst::m_e*std::sqrt( u2_sq ); p2x = PhysConst::m_e*u2x[j_2]; p2y = PhysConst::m_e*u2y[j_2]; p2z = PhysConst::m_e*u2z[j_2]; @@ -299,9 +303,9 @@ void DifferentialLuminosity2D::ComputeDiags (int step) Real const radicand = beta1_sq + beta2_sq - 2*beta1_dot_beta2 - beta1_sq*beta2_sq + beta1_dot_beta2*beta1_dot_beta2; Real const d2L_dE1_dE2 = PhysConst::c * std::sqrt( radicand ) * w1[j_1] * w2[j_2] / (dV * bin_size_1 * bin_size_2) * dt; // m^-2 eV^-2 - if (d2L_dE1_dE2 < 0){ - amrex::Print(0) << "NEGATIVE VALUE!!!!!!!!!!!!!!!!!!!!!!!!!\n"; - } + + //amrex::Print(0) << "DiffLumi 2D = " << d2L_dE1_dE2 << " bin size 1 = " << bin_size_1 << " bin size 2 = " << bin_size_2 << "\n"; + amrex::Real &data = d_table(bin_1, bin_2); amrex::HostDevice::Atomic::Add(&data, d2L_dE1_dE2); @@ -318,12 +322,12 @@ void DifferentialLuminosity2D::ComputeDiags (int step) if (m_intervals.contains(step+1)) { // Copy data from GPU memory - m_h_data_2D.copy(d_data_2D); + m_h_data_2D.copy(m_d_data_2D); // reduced sum over mpi ranks - const int size = static_cast (d_data_2D.size()); + const int size = static_cast (m_d_data_2D.size()); ParallelDescriptor::ReduceRealSum - (h_table_data.p, size, ParallelDescriptor::IOProcessorNumber()); + (m_h_data_2D.table().p, size, ParallelDescriptor::IOProcessorNumber()); } // Return for all that are not IO processor @@ -378,7 +382,7 @@ void DifferentialLuminosity2D::WriteToFile (int step) const // record components auto data = f_mesh[io::RecordComponent::SCALAR]; - data.setUnitSI(std::pow(PhysConst::q_e, -2)); + //data.setUnitSI(std::pow(PhysConst::q_e, -2)); // meta data f_mesh.setAxisLabels({"E2", "E1"}); // ordinate, abscissa From 46d51a220c991e8f9b86d10e095641f5045d27ea Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Mon, 13 Jan 2025 16:10:40 -0800 Subject: [PATCH 4/8] changed units and test --- Examples/Tests/diff_lumi_diag/inputs_base_3d | 17 ++--- .../ReducedDiags/DifferentialLuminosity2D.H | 2 +- .../ReducedDiags/DifferentialLuminosity2D.cpp | 67 +++++++------------ 3 files changed, 36 insertions(+), 50 deletions(-) diff --git a/Examples/Tests/diff_lumi_diag/inputs_base_3d b/Examples/Tests/diff_lumi_diag/inputs_base_3d index 4724fec3af0..5a70f60212f 100644 --- a/Examples/Tests/diff_lumi_diag/inputs_base_3d +++ b/Examples/Tests/diff_lumi_diag/inputs_base_3d @@ -28,6 +28,7 @@ my_constants.dt = sigmaz/clight/10. ################################# ####### GENERAL PARAMETERS ###### ################################# + stop_time = T amr.n_cell = nx ny nz amr.max_grid_size = 128 @@ -96,18 +97,18 @@ diag1.species = beam1 beam2 warpx.reduced_diags_names = DifferentialLuminosity_beam1_beam2 DifferentialLuminosity2d_beam1_beam2 DifferentialLuminosity_beam1_beam2.type = DifferentialLuminosity -DifferentialLuminosity_beam1_beam2.intervals = 5 +DifferentialLuminosity_beam1_beam2.intervals = 80 DifferentialLuminosity_beam1_beam2.species = beam1 beam2 DifferentialLuminosity_beam1_beam2.bin_number = 128 DifferentialLuminosity_beam1_beam2.bin_max = 2.1*beam_energy_eV -DifferentialLuminosity_beam1_beam2.bin_min = 1.9*beam_energy_eV +DifferentialLuminosity_beam1_beam2.bin_min = 0 DifferentialLuminosity2d_beam1_beam2.type = DifferentialLuminosity2D -DifferentialLuminosity2d_beam1_beam2.intervals = 5 +DifferentialLuminosity2d_beam1_beam2.intervals = 80 DifferentialLuminosity2d_beam1_beam2.species = beam1 beam2 DifferentialLuminosity2d_beam1_beam2.bin_number_1 = 128 -DifferentialLuminosity2d_beam1_beam2.bin_max_1 = 1.1*beam_energy_eV -DifferentialLuminosity2d_beam1_beam2.bin_min_1 = 0.9*beam_energy_eV -DifferentialLuminosity2d_beam1_beam2.bin_number_2 = 256 -DifferentialLuminosity2d_beam1_beam2.bin_max_2 = 1.1*beam_energy_eV -DifferentialLuminosity2d_beam1_beam2.bin_min_2 = 0.9*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_max_1 = 1.05*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_min_1 = 0 +DifferentialLuminosity2d_beam1_beam2.bin_number_2 = 128 +DifferentialLuminosity2d_beam1_beam2.bin_max_2 = 1.05*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_min_2 = 0 diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H index d3c776ee9f8..7ffefec324e 100644 --- a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H @@ -62,7 +62,7 @@ public: private: - /// output data // Array in which to accumulate the luminosity across timesteps + /// output table in which to accumulate the luminosity across timesteps amrex::TableData m_d_data_2D; }; diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp index b66e036934c..24691c5d32d 100644 --- a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp @@ -78,7 +78,7 @@ DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) // RZ coordinate is not supported #if (defined WARPX_DIM_RZ) WARPX_ABORT_WITH_MESSAGE( - "DifferentialLuminosity2D diagnostics do not work in RZ geometry."); + "DifferentialLuminosity2D diagnostics does not work in RZ geometry."); #endif // read colliding species names - must be 2 @@ -89,9 +89,6 @@ DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) m_beam_name.size() == 2u, "DifferentialLuminosity2D diagnostics must involve exactly two species"); - // should we add an assert that checks that the selected species actually exist? - // ref. ParticleHistrogram2D - pp_rd_name.query("openpmd_backend", m_openpmd_backend); pp_rd_name.query("file_min_digits", m_file_min_digits); // pick first available backend if default is chosen @@ -100,7 +97,7 @@ DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) } pp_rd_name.add("openpmd_backend", m_openpmd_backend); - // read bin parameters for first species + // read bin parameters for species 1 int bin_num_1 = 0; amrex::Real bin_max_1 = 0.0_rt, bin_min_1 = 0.0_rt; utils::parser::getWithParser(pp_rd_name, "bin_number_1", bin_num_1); @@ -111,7 +108,7 @@ DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) m_bin_min_1 = bin_min_1; m_bin_size_1 = (bin_max_1 - bin_min_1) / bin_num_1; - // read bin parameters for second species + // read bin parameters for species 2 int bin_num_2 = 0; amrex::Real bin_max_2 = 0.0_rt, bin_min_2 = 0.0_rt; utils::parser::getWithParser(pp_rd_name, "bin_number_2", bin_num_2); @@ -122,27 +119,25 @@ DifferentialLuminosity2D::DifferentialLuminosity2D (const std::string& rd_name) m_bin_min_2 = bin_min_2; m_bin_size_2 = (bin_max_2 - bin_min_2) / bin_num_2; - - // resize data array + // resize data array on the host Array tlo{0,0}; // lower bounds Array thi{m_bin_num_1-1, m_bin_num_2-1}; // inclusive upper bounds m_h_data_2D.resize(tlo, thi, The_Pinned_Arena()); - auto const& h_table_data = m_h_data_2D.table(); - - // Initialize data on the host + // initialize data on the host for (int i = tlo[0]; i <= thi[0]; ++i) { for (int j = tlo[1]; j <= thi[1]; ++j) { h_table_data(i,j) = 0.0_rt; } } + // resize data on the host m_d_data_2D.resize(tlo, thi); + // copy data from host to device m_d_data_2D.copy(m_h_data_2D); Gpu::streamSynchronize(); - -} +} // end constructor void DifferentialLuminosity2D::ComputeDiags (int step) { @@ -152,14 +147,15 @@ void DifferentialLuminosity2D::ComputeDiags (int step) WARPX_PROFILE("DifferentialLuminosity2D::ComputeDiags"); - auto d_table = m_d_data_2D.table(); - // Since this diagnostic *accumulates* the luminosity in the - // array d_data, we add contributions at *each timestep*, but + // table m_d_data_2D, we add contributions at *each timestep*, but // we only write the data to file at intervals specified by the user. const Real c_sq = PhysConst::c*PhysConst::c; const Real c_over_qe = PhysConst::c/PhysConst::q_e; + // output table data + auto d_table = m_d_data_2D.table(); + // get a reference to WarpX instance auto& warpx = WarpX::GetInstance(); const Real dt = warpx.getdt(0); @@ -202,18 +198,19 @@ void DifferentialLuminosity2D::ComputeDiags (int step) ParticleBins bins_1 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_1 ); ParticleBins bins_2 = ParticleUtils::findParticlesInEachCell( lev, mfi, ptile_2 ); - // Species + // species 1 const auto soa_1 = ptile_1.getParticleTileData(); index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr(); index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr(); - // Particle data in the tile/box + // extract particle data of species 1 in the current tile/box amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; - amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; - amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; // v*gamma=p/m + amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; // u=v*gamma=p/m + amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; bool const species1_is_photon = species_1.AmIA(); + // same for species 2 const auto soa_2 = ptile_2.getParticleTileData(); index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr(); index_type const* AMREX_RESTRICT cell_offsets_2 = bins_2.offsetsPtr(); @@ -224,7 +221,7 @@ void DifferentialLuminosity2D::ComputeDiags (int step) amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; bool const species2_is_photon = species_2.AmIA(); - // Extract low-level data + // Extract low-level (cell-level) data auto const n_cells = static_cast(bins_1.numBins()); // Loop over cells @@ -279,15 +276,14 @@ void DifferentialLuminosity2D::ComputeDiags (int step) Real const E_1 = p1t * c_over_qe; // eV Real const E_2 = p2t * c_over_qe; // eV - // determine particle 1 bin + // determine energy bin of particle 1 int const bin_1 = int(Math::floor((E_1-bin_min_1)/bin_size_1)); if ( bin_1<0 || bin_1>=num_bins_1 ) { continue; } // discard if out-of-range - // determine particle 2 bin + // determine energy bin of particle 2 int const bin_2 = int(Math::floor((E_2-bin_min_2)/bin_size_2)); if ( bin_2<0 || bin_2>=num_bins_2 ) { continue; } // discard if out-of-range - Real const inv_p1t = 1.0_rt/p1t; Real const inv_p2t = 1.0_rt/p2t; @@ -299,13 +295,10 @@ void DifferentialLuminosity2D::ComputeDiags (int step) // (v1 - v2)^2 = v1^2 + v2^2 - 2 v1.v2 // and (v1 x v2)^2 = v1^2 v2^2 - (v1.v2)^2 // we also use beta=v/c instead of v - Real const radicand = beta1_sq + beta2_sq - 2*beta1_dot_beta2 - beta1_sq*beta2_sq + beta1_dot_beta2*beta1_dot_beta2; Real const d2L_dE1_dE2 = PhysConst::c * std::sqrt( radicand ) * w1[j_1] * w2[j_2] / (dV * bin_size_1 * bin_size_2) * dt; // m^-2 eV^-2 - //amrex::Print(0) << "DiffLumi 2D = " << d2L_dE1_dE2 << " bin size 1 = " << bin_size_1 << " bin size 2 = " << bin_size_2 << "\n"; - amrex::Real &data = d_table(bin_1, bin_2); amrex::HostDevice::Atomic::Add(&data, d2L_dE1_dE2); @@ -315,7 +308,6 @@ void DifferentialLuminosity2D::ComputeDiags (int step) } // boxes } // levels - // Only write to file at intervals specified by the user. // At these intervals, the data needs to ready on the CPU, // so we copy it from the GPU to the CPU and reduce across MPI ranks. @@ -336,9 +328,6 @@ void DifferentialLuminosity2D::ComputeDiags (int step) #endif // not RZ } // end void DifferentialLuminosity2D::ComputeDiags - - - void DifferentialLuminosity2D::WriteToFile (int step) const { // Judge if the diags should be done at this step @@ -368,28 +357,26 @@ void DifferentialLuminosity2D::WriteToFile (int step) const io::Access::CREATE); auto i = series.iterations[step + 1]; // record - auto f_mesh = i.meshes["data"]; - //f_mesh.setAttribute("energy_1", "E1"); - //f_mesh.setAttribute("energy_2", "E2"); + auto f_mesh = i.meshes["d2L_dE1_dE2"]; f_mesh.setUnitDimension({ {io::UnitDimension::L, -6}, {io::UnitDimension::M, -2}, {io::UnitDimension::T, 4} }); - // add units! - // record components auto data = f_mesh[io::RecordComponent::SCALAR]; - //data.setUnitSI(std::pow(PhysConst::q_e, -2)); - // meta data f_mesh.setAxisLabels({"E2", "E1"}); // ordinate, abscissa std::vector< double > const& gridGlobalOffset = {m_bin_min_2, m_bin_min_1}; f_mesh.setGridGlobalOffset(gridGlobalOffset); f_mesh.setGridSpacing({m_bin_size_2, m_bin_size_1}); + // to save the data in SI add the following: + // f_mesh.setGridUnitSI(PhysConst::q_e); + // data.setUnitSI(std::pow(PhysConst::q_e, -2)); + data.setPosition({0.5, 0.5}); auto dataset = io::Dataset( @@ -397,8 +384,6 @@ void DifferentialLuminosity2D::WriteToFile (int step) const {static_cast(m_bin_num_2), static_cast(m_bin_num_1)}); data.resetDataset(dataset); - // UNIT DIMENSION IS NOT SET ON THE VALUES - // Get time at level 0 auto & warpx = WarpX::GetInstance(); auto const time = warpx.gett_new(0); @@ -415,6 +400,6 @@ void DifferentialLuminosity2D::WriteToFile (int step) const series.close(); #else amrex::ignore_unused(step); - WARPX_ABORT_WITH_MESSAGE("ParticleHistogram2D: Needs openPMD-api compiled into WarpX, but was not found!"); + WARPX_ABORT_WITH_MESSAGE("DifferentialLuminosity2D: Needs openPMD-api compiled into WarpX, but was not found!"); #endif } From c6564576fe0c24c224d328510ea97f65e1845a27 Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Mon, 13 Jan 2025 17:52:17 -0800 Subject: [PATCH 5/8] add docs --- Docs/source/usage/parameters.rst | 43 ++++++++++++++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 31f0e06ab5b..b2ba7d11202 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -3618,6 +3618,49 @@ This shifts analysis from post-processing to runtime calculation of reduction op * ``.bin_min`` (`float`, in eV) The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed. + * ``DifferentialLuminosity2D`` + This type computes the two-dimensional differential luminosity between two species, defined as: + + .. math:: + + \frac{d^2\mathcal{L}}{dE_1 dE_2}(E_1, E_2, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\; + \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(E_1 - E_1(\boldsymbol{p}_1)) \delta(E_2 - E_2(\boldsymbol{p}_2)) + + where :math:`f_i` is the distribution function of species :math:`i`, + :math:`\boldsymbol{p}_i` and :math:`E_i (\boldsymbol{p}_i) = \sqrt{m_1^2c^4 + c^2 |\boldsymbol{p}_i|^2}` + are, respectively, the momentum and the energy of a particle of the :math:`i`-th species. + The 2D differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-2}`. + + * ``.species`` (`list of two strings`) + The names of the two species for which the differential luminosity is computed. + + * ``.bin_number_1`` (`int` > 0) + The number of bins in energy :math:`E_1` + + * ``.bin_max_1`` (`float`, in eV) + The minimum value of :math:`E_1` for which the 2D differential luminosity is computed. + + * ``.bin_min_1`` (`float`, in eV) + The maximum value of :math:`E_2` for which the 2D differential luminosity is compute + + * ``.bin_number_2`` (`int` > 0) + The number of bins in energy :math:`E_2` + + * ``.bin_max_2`` (`float`, in eV) + The minimum value of :math:`E_2` for which the 2D differential luminosity is computed. + + * ``.bin_min_2`` (`float`, in eV) + The minimum value of :math:`E_2` for which the 2D differential luminosity is computed. + + * ``.file_min_digits`` (`int`) optional (default `6`) + The minimum number of digits used for the iteration number appended to the diagnostic file names. + + The output is a ```` folder containing a set of openPMD files. + The values of the diagnostic are stored in a record labeled `d2L_dE1_dE2`. + An example input file and a loading python script of + using the DifferentialLuminosity2D reduced diagnostics + are given in ``Examples/Tests/diff_lumi_diag/``. + * ``Timestep`` This type outputs the simulation's physical timestep (in seconds) at each mesh refinement level. From bfc77fb122ae8030387086a1f44a9d66c245d296 Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Mon, 13 Jan 2025 17:52:55 -0800 Subject: [PATCH 6/8] update test --- Examples/Tests/diff_lumi_diag/analysis.py | 48 ++++++++++++++------ Examples/Tests/diff_lumi_diag/inputs_base_3d | 4 +- 2 files changed, 37 insertions(+), 15 deletions(-) diff --git a/Examples/Tests/diff_lumi_diag/analysis.py b/Examples/Tests/diff_lumi_diag/analysis.py index cadb21023ab..1a810638d69 100755 --- a/Examples/Tests/diff_lumi_diag/analysis.py +++ b/Examples/Tests/diff_lumi_diag/analysis.py @@ -5,15 +5,20 @@ # In that case, the differential luminosity can be calculated analytically. import os +import re import numpy as np -from read_raw_data import read_reduced_diags_histogram +from openpmd_viewer import OpenPMDTimeSeries -# Extract the differential luminosity from the file -_, _, E_bin, bin_data = read_reduced_diags_histogram( - "./diags/reducedfiles/DifferentialLuminosity_beam1_beam2.txt" -) -dL_dE_sim = bin_data[-1] # Differential luminosity at the end of the simulation +# Extract the 1D differential luminosity from the file +filename = "./diags/reducedfiles/DifferentialLuminosity_beam1_beam2.txt" +with open(filename) as f: + # First line: header, contains the energies + line = f.readline() + E_bin = np.array(list(map(float, re.findall("=(.*?)\(", line)))) +data = np.loadtxt(filename) +dE_bin = E_bin[1] - E_bin[0] +dL_dE_sim = data[-1, 2:] # Differential luminosity at the end of the simulation # Beam parameters N = 1.2e10 @@ -33,21 +38,38 @@ * np.exp(-((E_bin - 2 * E_beam) ** 2) / (2 * sigma_E**2)) ) +# Extract the 2D differential luminosity from the file +series = OpenPMDTimeSeries("./diags/reducedfiles/DifferentialLuminosity2d_beam1_beam2/") +d2L_dE1_dE2, info = series.get_field("d2L_dE1_dE2", iteration=80) +E1, E2 = info.E1, info.E2 +dE1, dE2 = info.dE1, info.dE2 + # Extract test name from path test_name = os.path.split(os.getcwd())[1] print("test_name", test_name) # Pick tolerance if "leptons" in test_name: - tol = 1e-2 + tol1 = 1e-2 + tol2 = 1e-5 elif "photons" in test_name: # In the photons case, the particles are # initialized from a density distribution ; # tolerance is larger due to lower particle statistics - tol = 6e-2 + tol1 = 6e-2 + tol2 = 1e-5 + +# Check that the 1D diagnostic and analytical result match +error1 = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max() +print("Relative error: ", error1) +print("Tolerance: ", tol1) + +# Check that the 2D and 1D diagnostics match +error2 = abs(np.sum(d2L_dE1_dE2) * dE2 * dE1 - np.sum(dL_dE_sim) * dE_bin) / abs( + np.sum(dL_dE_sim) * dE_bin +) +print("Relative error: ", error2) +print("Tolerance: ", tol2) -# Check that the simulation result and analytical result match -error = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max() -print("Relative error: ", error) -print("Tolerance: ", tol) -assert error < tol +assert error1 < tol1 +assert error2 < tol2 diff --git a/Examples/Tests/diff_lumi_diag/inputs_base_3d b/Examples/Tests/diff_lumi_diag/inputs_base_3d index 5a70f60212f..0c65850e82b 100644 --- a/Examples/Tests/diff_lumi_diag/inputs_base_3d +++ b/Examples/Tests/diff_lumi_diag/inputs_base_3d @@ -107,8 +107,8 @@ DifferentialLuminosity2d_beam1_beam2.type = DifferentialLuminosity2D DifferentialLuminosity2d_beam1_beam2.intervals = 80 DifferentialLuminosity2d_beam1_beam2.species = beam1 beam2 DifferentialLuminosity2d_beam1_beam2.bin_number_1 = 128 -DifferentialLuminosity2d_beam1_beam2.bin_max_1 = 1.05*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_max_1 = 1.45*beam_energy_eV DifferentialLuminosity2d_beam1_beam2.bin_min_1 = 0 DifferentialLuminosity2d_beam1_beam2.bin_number_2 = 128 -DifferentialLuminosity2d_beam1_beam2.bin_max_2 = 1.05*beam_energy_eV +DifferentialLuminosity2d_beam1_beam2.bin_max_2 = 1.45*beam_energy_eV DifferentialLuminosity2d_beam1_beam2.bin_min_2 = 0 From 95aedcfe390e3242e85377bfc064d39af63572b4 Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Mon, 13 Jan 2025 19:38:30 -0800 Subject: [PATCH 7/8] cleanup docs --- Docs/source/usage/parameters.rst | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index b2ba7d11202..98192294a46 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -3623,10 +3623,13 @@ This shifts analysis from post-processing to runtime calculation of reduction op .. math:: - \frac{d^2\mathcal{L}}{dE_1 dE_2}(E_1, E_2, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\; - \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(E_1 - E_1(\boldsymbol{p}_1)) \delta(E_2 - E_2(\boldsymbol{p}_2)) + \frac{d^2\mathcal{L}}{dE_1 dE_2}(E_1, E_2, t) = \int_0^t dt'\int d\boldsymbol{x}\, \int d\boldsymbol{p}_1 \int d\boldsymbol{p}_2\; + \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ + f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(E_1 - E_1(\boldsymbol{p}_1)) \delta(E_2 - E_2(\boldsymbol{p}_2)) - where :math:`f_i` is the distribution function of species :math:`i`, + where :math:`f_i` is the distribution function of species :math:`i` + (normalized such that :math:`\int \int f(\boldsymbol{x} \boldsymbol{p}, t )d\boldsymbol{x} d\boldsymbol{p} = N` + is the number of particles in species :math:`i` at time :math:`t`), :math:`\boldsymbol{p}_i` and :math:`E_i (\boldsymbol{p}_i) = \sqrt{m_1^2c^4 + c^2 |\boldsymbol{p}_i|^2}` are, respectively, the momentum and the energy of a particle of the :math:`i`-th species. The 2D differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-2}`. From a5be100d5bbb03c9454b12df42ac17ae58fc8b24 Mon Sep 17 00:00:00 2001 From: Arianna Formenti Date: Tue, 14 Jan 2025 09:40:05 -0800 Subject: [PATCH 8/8] clean up and adjust tolerances --- Examples/Tests/diff_lumi_diag/analysis.py | 9 ++++----- .../ReducedDiags/DifferentialLuminosity2D.cpp | 8 ++------ 2 files changed, 6 insertions(+), 11 deletions(-) diff --git a/Examples/Tests/diff_lumi_diag/analysis.py b/Examples/Tests/diff_lumi_diag/analysis.py index 1a810638d69..9cb5bfdc0cf 100755 --- a/Examples/Tests/diff_lumi_diag/analysis.py +++ b/Examples/Tests/diff_lumi_diag/analysis.py @@ -41,7 +41,6 @@ # Extract the 2D differential luminosity from the file series = OpenPMDTimeSeries("./diags/reducedfiles/DifferentialLuminosity2d_beam1_beam2/") d2L_dE1_dE2, info = series.get_field("d2L_dE1_dE2", iteration=80) -E1, E2 = info.E1, info.E2 dE1, dE2 = info.dE1, info.dE2 # Extract test name from path @@ -50,14 +49,14 @@ # Pick tolerance if "leptons" in test_name: - tol1 = 1e-2 - tol2 = 1e-5 + tol1 = 0.02 + tol2 = 0.003 elif "photons" in test_name: # In the photons case, the particles are # initialized from a density distribution ; # tolerance is larger due to lower particle statistics - tol1 = 6e-2 - tol2 = 1e-5 + tol1 = 0.02 + tol2 = 0.003 # Check that the 1D diagnostic and analytical result match error1 = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max() diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp index 24691c5d32d..b3968b9fb02 100644 --- a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp +++ b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.cpp @@ -357,7 +357,7 @@ void DifferentialLuminosity2D::WriteToFile (int step) const io::Access::CREATE); auto i = series.iterations[step + 1]; // record - auto f_mesh = i.meshes["d2L_dE1_dE2"]; + auto f_mesh = i.meshes["d2L_dE1_dE2"]; // m^-2 eV^-2 f_mesh.setUnitDimension({ {io::UnitDimension::L, -6}, {io::UnitDimension::M, -2}, @@ -368,15 +368,11 @@ void DifferentialLuminosity2D::WriteToFile (int step) const auto data = f_mesh[io::RecordComponent::SCALAR]; // meta data - f_mesh.setAxisLabels({"E2", "E1"}); // ordinate, abscissa + f_mesh.setAxisLabels({"E2", "E1"}); // eV, eV std::vector< double > const& gridGlobalOffset = {m_bin_min_2, m_bin_min_1}; f_mesh.setGridGlobalOffset(gridGlobalOffset); f_mesh.setGridSpacing({m_bin_size_2, m_bin_size_1}); - // to save the data in SI add the following: - // f_mesh.setGridUnitSI(PhysConst::q_e); - // data.setUnitSI(std::pow(PhysConst::q_e, -2)); - data.setPosition({0.5, 0.5}); auto dataset = io::Dataset(