From 7af2163731a3895132b4bd87c5a735c79b3db5a1 Mon Sep 17 00:00:00 2001 From: Severin Date: Thu, 1 Jun 2023 17:40:55 +0200 Subject: [PATCH 01/11] add plasma insitu diags --- src/Hipace.cpp | 2 + src/particles/plasma/MultiPlasma.H | 8 + src/particles/plasma/MultiPlasma.cpp | 22 +- .../plasma/PlasmaParticleContainer.H | 44 +++- .../plasma/PlasmaParticleContainer.cpp | 245 +++++++++++++++++- 5 files changed, 318 insertions(+), 3 deletions(-) diff --git a/src/Hipace.cpp b/src/Hipace.cpp index 7fbf8fce56..a2210e090b 100644 --- a/src/Hipace.cpp +++ b/src/Hipace.cpp @@ -459,6 +459,7 @@ Hipace::Evolve () } m_multi_beam.InSituWriteToFile(step, m_physical_time, m_3D_geom[0]); + m_multi_plasma.InSituWriteToFile(step, m_physical_time, m_3D_geom[0]); // printing and resetting predictor corrector loop diagnostics if (m_verbose>=2 && !m_explicit) amrex::AllPrint() << "Rank " << rank << @@ -485,6 +486,7 @@ Hipace::SolveOneSlice (int islice, const int islice_local, int step) amrex::ParallelContext::push(m_comm_xy); m_multi_beam.InSituComputeDiags(step, islice, islice_local); + m_multi_plasma.InSituComputeDiags(step, islice); // Get this laser slice from the 3D array m_multi_laser.Copy(islice, false); diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index 01403ed12b..cc4bf9ff51 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -143,6 +143,14 @@ public: */ void TagByLevel (const int current_N_level, amrex::Vector const& geom3D); + /** Compute reduced plasma diagnostics of current slice, store in member variable. + * \param[in] step time step of simulation + * \param[in] islice current slice, on which diags are computed. + * \param[in] islice_local local index of the slice + */ + void InSituComputeDiags (int step, int islice); + void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom); + int m_sort_bin_size {32}; /**< Tile size to sort plasma particles */ /** Number of binary collisions */ diff --git a/src/particles/plasma/MultiPlasma.cpp b/src/particles/plasma/MultiPlasma.cpp index b84a87f52f..259f1fafe6 100644 --- a/src/particles/plasma/MultiPlasma.cpp +++ b/src/particles/plasma/MultiPlasma.cpp @@ -50,7 +50,7 @@ MultiPlasma::InitData (amrex::Vector slice_ba, for (auto& plasma : m_all_plasmas) { // make it think there is only level 0 plasma.SetParGDB(slice_gm[0], slice_dm[0], slice_ba[0]); - plasma.InitData(); + plasma.InitData(gm[0]); if(plasma.m_can_ionize) { PlasmaParticleContainer* plasma_product = nullptr; @@ -196,3 +196,23 @@ MultiPlasma::TagByLevel (const int current_N_level, amrex::Vector m_density_func; /**< Density function for the plasma */ amrex::Real m_min_density {0.}; /**< minimal density at which particles are injected */ @@ -171,6 +184,35 @@ public: private: std::string m_name; /**< name of the species */ + int m_nslices; /**< number of z slices of the domain */ + /** Number of real plasma properties for in-situ per-slice reduced diagnostics. */ + int m_insitu_nrp {13}; + /** Number of int plasma properties for in-situ per-slice reduced diagnostics. */ + int m_insitu_nip {1}; + /** Per-slice real plasma properties: + * 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 + * sum(w), [x], [x^2], [y], [y^2], [ux], [ux^2], [uy], [uy^2], [x*ux], [y*uy], [ga], [ga^2] + * where [] means average over all particles within slice. + * Per-slice emittance: sqrt( abs( ([x^2]-[x]^2) * ([ux^2]-[ux]^2) - ([x*ux]-[x][ux])^2 ) ). + * Projected emittance: Same as above AFTER averaging all these quantities over slices. + * Energy spread: sqrt([ga^2]-[ga]^2), and same as above. + * Np: number of particles in this slice + */ + amrex::Vector m_insitu_rdata; + /** Per-slice int plasma properties: + * 0 + * Np + * Np: number of particles in this slice + */ + amrex::Vector m_insitu_idata; + /** Sum of all per-slice real plasma properties */ + amrex::Vector m_insitu_sum_rdata; + /** Sum of all per-slice int plasma properties */ + amrex::Vector m_insitu_sum_idata; + /** Prefix/path for the output files */ + std::string m_insitu_file_prefix = "diags/plasma_insitu"; + /** How often the insitu plasma diagnostics should be computed and written */ + int m_insitu_period {-1}; }; /** \brief Iterator over boxes in a particle container */ diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index 1cf6500126..9d10ece2c5 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -12,6 +12,10 @@ #include "utils/AtomicWeightTable.H" #include "utils/DeprecatedInput.H" #include "utils/GPUUtil.H" +#include "utils/InsituUtil.H" +#ifdef HIPACE_USE_OPENPMD +# include +#endif #include "particles/pusher/PlasmaParticleAdvance.H" #include "particles/pusher/BeamParticleAdvance.H" #include "particles/particles_utils/FieldGather.H" @@ -155,15 +159,32 @@ PlasmaParticleContainer::ReadParameters () {Hipace::m_depos_order_xy % 2, Hipace::m_depos_order_xy % 2}; queryWithParserAlt(pp, "reorder_idx_type", idx_array, pp_alt); m_reorder_idx_type = amrex::IntVect(idx_array[0], idx_array[1], 0); + queryWithParserAlt(pp, "insitu_period", m_insitu_period, pp_alt); + queryWithParserAlt(pp, "insitu_file_prefix", m_insitu_file_prefix, pp_alt); } void -PlasmaParticleContainer::InitData () +PlasmaParticleContainer::InitData (const amrex::Geometry& geom) { reserveData(); resizeData(); InitParticles(m_ppc, m_u_std, m_u_mean, m_radius, m_hollow_core_radius); + + if (m_insitu_period > 0) { + // TODO FIXME +// #ifdef HIPACE_USE_OPENPMD +// AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_insitu_file_prefix != +// Hipace::GetInstance().m_openpmd_writer.m_file_prefix, +// "Must choose a different insitu file prefix compared to the full diagnostics"); +// #endif + // Allocate memory for in-situ diagnostics + m_nslices = geom.Domain().length(2); + m_insitu_rdata.resize(m_nslices*m_insitu_nrp, 0.); + m_insitu_idata.resize(m_nslices*m_insitu_nip, 0); + m_insitu_sum_rdata.resize(m_insitu_nrp, 0.); + m_insitu_sum_idata.resize(m_insitu_nip, 0); + } } void @@ -410,3 +431,225 @@ IonizationModule (const int lev, }); } } + +bool +PlasmaParticleContainer::doInSitu (int step) +{ + if (m_insitu_period <= 0) return false; + return step % m_insitu_period == 0; +} + +void +PlasmaParticleContainer::InSituComputeDiags (int islice) +{ + HIPACE_PROFILE("PlasmaParticleContainer::InSituComputeDiags()"); + + using namespace amrex::literals; + + AMREX_ALWAYS_ASSERT(m_insitu_rdata.size()>0 && m_insitu_idata.size()>0 && + m_insitu_sum_rdata.size()>0 && m_insitu_sum_idata.size()>0); + + PhysConst const phys_const = get_phys_const(); + const amrex::Real clight_inv = 1.0_rt/phys_const.c; + const amrex::Real clightsq_inv = 1.0_rt/(phys_const.c*phys_const.c); + + // Loop over particle boxes + for (PlasmaParticleIterator pti(*this); pti.isValid(); ++pti) + { + // loading the data + const auto ptd = pti.GetParticleTile().getParticleTileData(); + + const int num_particles = pti.numParticles(); + // #ifdef AMREX_USE_OMP + // #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) + // #endif + // { + // amrex::Long const num_particles = pti.numParticles(); + // #ifdef AMREX_USE_OMP + // amrex::Long const idx_begin = (num_particles * omp_get_thread_num()) / omp_get_num_threads(); + // amrex::Long const idx_end = (num_particles * (omp_get_thread_num()+1)) / omp_get_num_threads(); + // #else + // amrex::Long constexpr idx_begin = 0; + // amrex::Long const idx_end = num_particles; + // #endif + // Tuple contains: + // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 + // sum(w), , , , , , , , , , , , , np + amrex::ReduceOps reduce_op; + amrex::ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + reduce_op.eval( + num_particles, reduce_data, + [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple + { + const int ip = i; //indices[cell_start+i]; + + const amrex::Real xp = ptd.pos(0, ip); + const amrex::Real yp = ptd.pos(1, ip); + const amrex::Real uxp = ptd.rdata(PlasmaIdx::ux )[ip] * clight_inv; // proper velocity to u + const amrex::Real uyp = ptd.rdata(PlasmaIdx::uy )[ip] * clight_inv; + const amrex::Real psip = ptd.rdata(PlasmaIdx::psi)[ip]; + const amrex::Real wp = ptd.rdata(PlasmaIdx::w )[ip]; + + if (ptd.id(ip) < 0) { + return{0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, + 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0._rt, 0}; + } + + // particle's Lorentz factor + const amrex::Real gamma = (1.0_rt + uxp*uxp*clightsq_inv + + uyp*uyp*clightsq_inv + psip*psip)/(2.0_rt*psip); + amrex::Real uzp = (gamma - psip); // the *c from uz cancels with the /c from the proper velocity conversion + + return {wp, + wp*xp, + wp*xp*xp, + wp*yp, + wp*yp*yp, + wp*uxp, + wp*uxp*uxp, + wp*uyp, + wp*uyp*uyp, + wp*uzp, + wp*uzp*uzp, + wp*gamma, + wp*gamma*gamma, + 1}; + }); + + ReduceTuple a = reduce_data.value(); + const amrex::Real sum_w0 = amrex::get< 0>(a); + const amrex::Real sum_w_inv = sum_w0::epsilon() ? 0._rt : 1._rt/sum_w0; + + m_insitu_rdata[islice ] = sum_w0; + m_insitu_rdata[islice+ 1*m_nslices] = amrex::get< 1>(a)*sum_w_inv; + m_insitu_rdata[islice+ 2*m_nslices] = amrex::get< 2>(a)*sum_w_inv; + m_insitu_rdata[islice+ 3*m_nslices] = amrex::get< 3>(a)*sum_w_inv; + m_insitu_rdata[islice+ 4*m_nslices] = amrex::get< 4>(a)*sum_w_inv; + m_insitu_rdata[islice+ 5*m_nslices] = amrex::get< 5>(a)*sum_w_inv; + m_insitu_rdata[islice+ 6*m_nslices] = amrex::get< 6>(a)*sum_w_inv; + m_insitu_rdata[islice+ 7*m_nslices] = amrex::get< 7>(a)*sum_w_inv; + m_insitu_rdata[islice+ 8*m_nslices] = amrex::get< 8>(a)*sum_w_inv; + m_insitu_rdata[islice+ 9*m_nslices] = amrex::get< 9>(a)*sum_w_inv; + m_insitu_rdata[islice+10*m_nslices] = amrex::get<10>(a)*sum_w_inv; + m_insitu_rdata[islice+11*m_nslices] = amrex::get<11>(a)*sum_w_inv; + m_insitu_rdata[islice+12*m_nslices] = amrex::get<12>(a)*sum_w_inv; + m_insitu_idata[islice ] = amrex::get<13>(a); + + m_insitu_sum_rdata[ 0] += sum_w0; + m_insitu_sum_rdata[ 1] += amrex::get< 1>(a); + m_insitu_sum_rdata[ 2] += amrex::get< 2>(a); + m_insitu_sum_rdata[ 3] += amrex::get< 3>(a); + m_insitu_sum_rdata[ 4] += amrex::get< 4>(a); + m_insitu_sum_rdata[ 5] += amrex::get< 5>(a); + m_insitu_sum_rdata[ 6] += amrex::get< 6>(a); + m_insitu_sum_rdata[ 7] += amrex::get< 7>(a); + m_insitu_sum_rdata[ 8] += amrex::get< 8>(a); + m_insitu_sum_rdata[ 9] += amrex::get< 9>(a); + m_insitu_sum_rdata[10] += amrex::get<10>(a); + m_insitu_sum_rdata[11] += amrex::get<11>(a); + m_insitu_sum_rdata[12] += amrex::get<12>(a); + m_insitu_sum_idata[ 0] += amrex::get<13>(a); + } +} + +void +PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom) +{ + HIPACE_PROFILE("PlasmaParticleContainer::InSituWriteToFile()"); + +#ifdef HIPACE_USE_OPENPMD + // create subdirectory + openPMD::auxiliary::create_directories(m_insitu_file_prefix); +#endif + + // zero pad the rank number; + std::string::size_type n_zeros = 4; + std::string rank_num = std::to_string(amrex::ParallelDescriptor::MyProc()); + std::string pad_rank_num = std::string(n_zeros-std::min(rank_num.size(), n_zeros),'0')+rank_num; + + // open file + std::ofstream ofs{m_insitu_file_prefix + "/reduced_" + m_name + "." + pad_rank_num + ".txt", + std::ofstream::out | std::ofstream::app | std::ofstream::binary}; + + const amrex::Real sum_w0 = m_insitu_sum_rdata[0]; + const std::size_t nslices = static_cast(m_nslices); + const amrex::Real normalized_density_factor = Hipace::m_normalized_units ? + geom.CellSizeArray().product() : 1; // dx * dy * dz in normalized units, 1 otherwise + + // specify the structure of the data later available in python + // avoid pointers to temporary objects as second argument, stack variables are ok + const amrex::Vector all_data{ + {"time" , &time}, + {"step" , &step}, + {"n_slices", &m_nslices}, + {"charge" , &m_charge}, + {"mass" , &m_mass}, + {"z_lo" , &geom.ProbLo()[2]}, + {"z_hi" , &geom.ProbHi()[2]}, + {"normalized_density_factor", &normalized_density_factor}, + {"[x]" , &m_insitu_rdata[1*nslices], nslices}, + {"[x^2]" , &m_insitu_rdata[2*nslices], nslices}, + {"[y]" , &m_insitu_rdata[3*nslices], nslices}, + {"[y^2]" , &m_insitu_rdata[4*nslices], nslices}, + {"[ux]" , &m_insitu_rdata[5*nslices], nslices}, + {"[ux^2]" , &m_insitu_rdata[6*nslices], nslices}, + {"[uy]" , &m_insitu_rdata[7*nslices], nslices}, + {"[uy^2]" , &m_insitu_rdata[8*nslices], nslices}, + {"[uz]" , &m_insitu_rdata[9*nslices], nslices}, + {"[uz^2]" , &m_insitu_rdata[10*nslices], nslices}, + {"[ga]" , &m_insitu_rdata[11*nslices], nslices}, + {"[ga^2]" , &m_insitu_rdata[12*nslices], nslices}, + {"sum(w)" , &m_insitu_rdata[0], nslices}, + {"Np" , &m_insitu_idata[0], nslices}, + {"average" , { + {"[x]" , &(m_insitu_sum_rdata[ 1] /= sum_w0)}, + {"[x^2]" , &(m_insitu_sum_rdata[ 2] /= sum_w0)}, + {"[y]" , &(m_insitu_sum_rdata[ 3] /= sum_w0)}, + {"[y^2]" , &(m_insitu_sum_rdata[ 4] /= sum_w0)}, + {"[ux]" , &(m_insitu_sum_rdata[ 5] /= sum_w0)}, + {"[ux^2]", &(m_insitu_sum_rdata[ 6] /= sum_w0)}, + {"[uy]" , &(m_insitu_sum_rdata[ 7] /= sum_w0)}, + {"[uy^2]", &(m_insitu_sum_rdata[ 8] /= sum_w0)}, + {"[uz]", &(m_insitu_sum_rdata[ 9] /= sum_w0)}, + {"[uz^2]", &(m_insitu_sum_rdata[10] /= sum_w0)}, + {"[ga]" , &(m_insitu_sum_rdata[11] /= sum_w0)}, + {"[ga^2]", &(m_insitu_sum_rdata[12] /= sum_w0)} + }}, + {"total" , { + {"sum(w)", &m_insitu_sum_rdata[0]}, + {"Np" , &m_insitu_sum_idata[0]} + }} + }; + + if (ofs.tellp() == 0) { + // write JSON header containing a NumPy structured datatype + insitu_utils::write_header(all_data, ofs); + } + + // write binary data according to datatype in header + insitu_utils::write_data(all_data, ofs); + + // close file + ofs.close(); + // assert no file errors +#ifdef HIPACE_USE_OPENPMD + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ofs, "Error while writing insitu plasma diagnostics"); +#else + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ofs, "Error while writing insitu plasma diagnostics. " + "Maybe the specified subdirectory does not exist"); +#endif + + // reset arrays for insitu data + for (auto& x : m_insitu_rdata) x = 0.; + for (auto& x : m_insitu_idata) x = 0; + for (auto& x : m_insitu_sum_rdata) x = 0.; + for (auto& x : m_insitu_sum_idata) x = 0; +} From 7128f66fbcfe3623b3932b5fb687c5b2fe783d56 Mon Sep 17 00:00:00 2001 From: Severin Date: Fri, 2 Jun 2023 18:14:10 +0200 Subject: [PATCH 02/11] add preliminary doc --- docs/source/run/parameters.rst | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index e87fb2fee9..9985713a81 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -357,6 +357,40 @@ When both are specified, the per-species value is used. the specific ionization energy of each state. Options are: ``electron``, ``positron``, ``H``, ``D``, ``T``, ``He``, ``Li``, ``Be``, ``B``, …. +* `` or plasmas.insitu_period`` (`int`) optional (default ``-1``) + Period of in-situ diagnostics, for computing the main per-slice plasma quantities + (energy, temperature, etc.). This works in the same way as the beam insitu diagnostics. + For this the following quantities are calculated per slice and stored: + ``sum(w), [x], [x^2], [y], [y^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], [ga], [ga^2], np`` + where "[]" stands for averaging over all particles in the current slice, + "w" stands for weight, "ux" is the momentum in the x direction, "ga" is the Lorentz factor. + Averages and totals over all slices are also provided for convenience under the + respective ``average`` and ``total`` subcategories. + + Additionally, some metadata is also available: + ``time, step, n_slices, charge, mass, z_lo, z_hi, normalized_density_factor``. + ``time`` and ``step`` refers to the physical time of the simulation and step number of the + current time step. + ``n_slices`` equal to the number of slices in the zeta direction. + ``charge`` and ``mass`` relate to a single plasma particle and are for example equal to the + electron charge and mass. + ``z_lo`` and ``z_hi`` are the lower and upper bounds of the z-axis of the simulation domain + specified in the input file and can be used to generate a z/zeta-axis for plotting. + ``normalized_density_factor`` is equal to ``dx * dy * dz`` in normalized units and 1 in + SI units. It can be used to convert ``sum(w)``, which specifies the beam density in normalized + units and beam weight an SI units, to the beam weight in both unit systems. + + The data is written to a file at ``/reduced_..txt``. + The in-situ diagnostics file format consists of a header part in ASCII containing a JSON object. + When this is parsed into Python it can be converted to a NumPy structured datatype. + The rest of the file, following immediately after the closing }, is in binary format and + contains all of the in-situ diagnostic along with some meta data. This part can be read using the + structured datatype of the first section. + Use ``hipace/tools/read_insitu_diagnostics.py`` to read the files using this format. + +* `` or plasmas.insitu_file_prefix`` (`string`) optional (default ``"plasma_diags/insitu"``) + Path of the plasma in-situ output. + * ``.can_ionize`` (`bool`) optional (default `0`) Whether this plasma can ionize. Can also be set to 1 by specifying ``.ionization_product``. From f73bd4f7f029117cfe7350ae2f010555e5952401 Mon Sep 17 00:00:00 2001 From: Severin Date: Fri, 2 Jun 2023 18:18:25 +0200 Subject: [PATCH 03/11] fix doxygen --- src/particles/plasma/MultiPlasma.H | 1 - 1 file changed, 1 deletion(-) diff --git a/src/particles/plasma/MultiPlasma.H b/src/particles/plasma/MultiPlasma.H index cc4bf9ff51..02519625f4 100644 --- a/src/particles/plasma/MultiPlasma.H +++ b/src/particles/plasma/MultiPlasma.H @@ -146,7 +146,6 @@ public: /** Compute reduced plasma diagnostics of current slice, store in member variable. * \param[in] step time step of simulation * \param[in] islice current slice, on which diags are computed. - * \param[in] islice_local local index of the slice */ void InSituComputeDiags (int step, int islice); void InSituWriteToFile (int step, amrex::Real time, const amrex::Geometry& geom); From 9fd0912196cbe611063899862f8a3c056effedc3 Mon Sep 17 00:00:00 2001 From: Severin Date: Fri, 2 Jun 2023 22:43:51 +0200 Subject: [PATCH 04/11] improve doc, add assert, add functions to tools --- docs/source/run/parameters.rst | 123 ++++++++---------- src/particles/beam/BeamParticleContainer.cpp | 2 +- .../plasma/PlasmaParticleContainer.cpp | 15 ++- tools/read_insitu_diagnostics.py | 36 +++++ 4 files changed, 99 insertions(+), 77 deletions(-) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index a0331b03f6..686dab9ff6 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -357,40 +357,6 @@ When both are specified, the per-species value is used. the specific ionization energy of each state. Options are: ``electron``, ``positron``, ``H``, ``D``, ``T``, ``He``, ``Li``, ``Be``, ``B``, …. -* `` or plasmas.insitu_period`` (`int`) optional (default ``-1``) - Period of in-situ diagnostics, for computing the main per-slice plasma quantities - (energy, temperature, etc.). This works in the same way as the beam insitu diagnostics. - For this the following quantities are calculated per slice and stored: - ``sum(w), [x], [x^2], [y], [y^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], [ga], [ga^2], np`` - where "[]" stands for averaging over all particles in the current slice, - "w" stands for weight, "ux" is the momentum in the x direction, "ga" is the Lorentz factor. - Averages and totals over all slices are also provided for convenience under the - respective ``average`` and ``total`` subcategories. - - Additionally, some metadata is also available: - ``time, step, n_slices, charge, mass, z_lo, z_hi, normalized_density_factor``. - ``time`` and ``step`` refers to the physical time of the simulation and step number of the - current time step. - ``n_slices`` equal to the number of slices in the zeta direction. - ``charge`` and ``mass`` relate to a single plasma particle and are for example equal to the - electron charge and mass. - ``z_lo`` and ``z_hi`` are the lower and upper bounds of the z-axis of the simulation domain - specified in the input file and can be used to generate a z/zeta-axis for plotting. - ``normalized_density_factor`` is equal to ``dx * dy * dz`` in normalized units and 1 in - SI units. It can be used to convert ``sum(w)``, which specifies the beam density in normalized - units and beam weight an SI units, to the beam weight in both unit systems. - - The data is written to a file at ``/reduced_..txt``. - The in-situ diagnostics file format consists of a header part in ASCII containing a JSON object. - When this is parsed into Python it can be converted to a NumPy structured datatype. - The rest of the file, following immediately after the closing }, is in binary format and - contains all of the in-situ diagnostic along with some meta data. This part can be read using the - structured datatype of the first section. - Use ``hipace/tools/read_insitu_diagnostics.py`` to read the files using this format. - -* `` or plasmas.insitu_file_prefix`` (`string`) optional (default ``"plasma_diags/insitu"``) - Path of the plasma in-situ output. - * ``.can_ionize`` (`bool`) optional (default `0`) Whether this plasma can ionize. Can also be set to 1 by specifying ``.ionization_product``. @@ -525,40 +491,6 @@ which are valid only for certain beam types, are introduced further below under ``n_subcycles`` times with a time step of `dt/n_subcycles`. This can be used to improve accuracy in highly non-linear focusing fields. -* `` or beams.insitu_period`` (`int`) optional (default ``-1``) - Period of in-situ diagnostics, for computing the main per-slice beam quantities for the main - beam parameters (width, energy spread, emittance, etc.). - For this the following quantities are calculated per slice and stored: - ``sum(w), [x], [x^2], [y], [y^2], [z], [z^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], [x*ux], [y*uy], [z*uz], [ga], [ga^2], np`` - where "[]" stands for averaging over all particles in the current slice, - "w" stands for weight, "ux" is the momentum in the x direction, "ga" is the Lorentz factor. - Averages and totals over all slices are also provided for convenience under the - respective ``average`` and ``total`` subcategories. - - Additionally, some metadata is also available: - ``time, step, n_slices, charge, mass, z_lo, z_hi, normalized_density_factor``. - ``time`` and ``step`` refers to the physical time of the simulation and step number of the - current timestep. - ``n_slices`` equal to the number of slices in the zeta direction. - ``charge`` and ``mass`` relate to a single beam particle and are for example equal to the - electron charge and mass. - ``z_lo`` and ``z_hi`` are the lower and upper bounds of the z-axis of the simulation domain - specified in the input file and can be used to generate a z/zeta-axis for plotting. - ``normalized_density_factor`` is equal to ``dx * dy * dz`` in normalized units and 1 in - SI units. It can be used to convert ``sum(w)``, which specifies the beam density in normalized - units and beam weight an SI units, to the beam weight in both unit systems. - - The data is written to a file at ``/reduced_..txt``. - The in-situ diagnostics file format consists of a header part in ASCII containing a JSON object. - When this is parsed into Python it can be converted to a NumPy structured datatype. - The rest of the file, following immediately after the closing }, is in binary format and - contains all of the in-situ diagnostic along with some meta data. This part can be read using the - structured datatype of the first section. - Use ``hipace/tools/read_insitu_diagnostics.py`` to read the files using this format. - -* `` or beams.insitu_file_prefix`` (`string`) optional (default ``"diags/insitu"``) - Path of the in-situ output. - * ``.do_salame`` (`bool`) optional (default `0`) Whether to use the SALAME algorithm [S. Diederichs et al., Phys. Rev. Accel. Beams 23, 121301 (2020)] to automatically flatten the accelerating field in the first time step. If turned on, the per-slice beam weight in the first time-step is adjusted such that the Ez field will be uniform in the beam. @@ -720,8 +652,11 @@ Parameters starting with ``lasers.`` apply to all laser pulses, parameters start Diagnostic parameters --------------------- +There are different types of diagnostics in HiPACE++. The standard diagnostics are compliant with the openPMD standard. The +in-situ diagnostics allow for fast analysis of large beams or the plasma particles. + * ``diagnostic.output_period`` (`integer`) optional (default `0`) - Output period for all diagnostics. Field or beam specific diagnostics can overwrite this parameter. + Output period for standard beam and field diagnostics. Field or beam specific diagnostics can overwrite this parameter. No output is given for ``diagnostic.output_period = 0``. Beam diagnostics @@ -783,3 +718,53 @@ Field diagnostics * `` or diagnostic.patch_hi`` (3 `float`) optional (default `infinity infinity infinity`) Upper limit for the diagnostic grid. + +In-situ diagnostics +^^^^^^^^^^^^^^^^^^^ + +Besides the standard diagnostics, fast in-situ diagnostics are available. They are most useful when beams with a large number of particles is used, because then the important moments can be calculated by GPU and significantly reduce the post-analysis of the simulation. +For particle beams, they can be used to calculate the main characterizing beam parameters (width, energy spread, emittance, etc.). Additionally, the plasma particle properties (e.g, the temperature) can be calculated. + +The in-situ quantities are calculated per slice and together with the weights per slice, all important quantities can be calculated. + +For particle beams, the following quantities are calculated per slice and stored: +``sum(w), [x], [x^2], [y], [y^2], [z], [z^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], [x*ux], [y*uy], [z*uz], [ga], [ga^2], np``. +For plasma particles, the following quantities are calculated per slice and stored: +``sum(w), [x], [x^2], [y], [y^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], [ga], [ga^2], np``. +Thereby, "[]" stands for averaging over all particles in the current slice, +"w" stands for weight, "ux" is the normalized momentum in the x direction, "ga" is the Lorentz factor. +Averages and totals over all slices are also provided for convenience under the +respective ``average`` and ``total`` subcategories. + +Additionally, some metadata is also available: +``time, step, n_slices, charge, mass, z_lo, z_hi, normalized_density_factor``. +``time`` and ``step`` refers to the physical time of the simulation and step number of the +current timestep. +``n_slices`` equal to the number of slices in the zeta direction. +``charge`` and ``mass`` relate to a single particle and are for example equal to the +electron charge and mass. +``z_lo`` and ``z_hi`` are the lower and upper bounds of the z-axis of the simulation domain +specified in the input file and can be used to generate a z/zeta-axis for plotting. +``normalized_density_factor`` is equal to ``dx * dy * dz`` in normalized units and 1 in +SI units. It can be used to convert ``sum(w)``, which specifies the beam density in normalized +units and beam weight an SI units, to the beam weight in both unit systems. + +The data is written to a file at ``/reduced_..txt``. +The in-situ diagnostics file format consists of a header part in ASCII containing a JSON object. +When this is parsed into Python it can be converted to a NumPy structured datatype. +The rest of the file, following immediately after the closing }, is in binary format and +contains all of the in-situ diagnostic along with some meta data. This part can be read using the +structured datatype of the first section. +Use ``hipace/tools/read_insitu_diagnostics.py`` to read the files using this format. Functions to calculate the most useful properties are also provided in that file. + +* `` or beams.insitu_period`` (`int`) optional (default ``-1``) + Period of the beam in-situ diagnostics. `-1` means no in-situ diagnostics. + +* `` or beams.insitu_file_prefix`` (`string`) optional (default ``"diags/insitu"``) + Path of the beam in-situ output. Must not be the same as `hipace.file_prefix`. + +* `` or plasmas.insitu_period`` (`int`) optional (default ``-1``) + Period of the plasma in-situ diagnostics. `-1` means no in-situ diagnostics. + +* `` or plasmas.insitu_file_prefix`` (`string`) optional (default ``"plasma_diags/insitu"``) + Path of the plasma in-situ output. Must not be the same as `hipace.file_prefix`. diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp index d4fa1e86d2..71c472f81e 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -185,7 +185,7 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom) #ifdef HIPACE_USE_OPENPMD AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_insitu_file_prefix != Hipace::GetInstance().m_openpmd_writer.m_file_prefix, - "Must choose a different insitu file prefix compared to the full diagnostics"); + "Must choose a different beam insitu file prefix compared to the full diagnostics"); #endif // Allocate memory for in-situ diagnostics m_nslices = geom.Domain().length(2); diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index a6348e247b..ef575d7a01 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -170,12 +170,11 @@ PlasmaParticleContainer::InitData (const amrex::Geometry& geom) InitParticles(m_ppc, m_u_std, m_u_mean, m_radius, m_hollow_core_radius); if (m_insitu_period > 0) { - // TODO FIXME -// #ifdef HIPACE_USE_OPENPMD -// AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_insitu_file_prefix != -// Hipace::GetInstance().m_openpmd_writer.m_file_prefix, -// "Must choose a different insitu file prefix compared to the full diagnostics"); -// #endif +#ifdef HIPACE_USE_OPENPMD + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_insitu_file_prefix != + Hipace::GetInstance().m_openpmd_writer.m_file_prefix, + "Must choose a different plasma insitu file prefix compared to the full diagnostics"); +#endif // Allocate memory for in-situ diagnostics m_nslices = geom.Domain().length(2); m_insitu_rdata.resize(m_nslices*m_insitu_nrp, 0.); @@ -404,7 +403,7 @@ IonizationModule (const int lev, if(p_ion_mask[ip] != 0) { const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); - const long pidx = pid + old_size; + const long pidx = pid + old_size;beam // Copy ion data to new electron int_arrdata_elec[PlasmaIdx::id ][pidx] = pid_start + pid; @@ -588,6 +587,7 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am const std::size_t nslices = static_cast(m_nslices); const amrex::Real normalized_density_factor = Hipace::m_normalized_units ? geom.CellSizeArray().product() : 1; // dx * dy * dz in normalized units, 1 otherwise + const int is_normalized_units = Hipace::m_normalized_units; // specify the structure of the data later available in python // avoid pointers to temporary objects as second argument, stack variables are ok @@ -600,6 +600,7 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am {"z_lo" , &geom.ProbLo()[2]}, {"z_hi" , &geom.ProbHi()[2]}, {"normalized_density_factor", &normalized_density_factor}, + {"is_normalized_units", &is_normalized_units}, {"[x]" , &m_insitu_rdata[1*nslices], nslices}, {"[x^2]" , &m_insitu_rdata[2*nslices], nslices}, {"[y]" , &m_insitu_rdata[3*nslices], nslices}, diff --git a/tools/read_insitu_diagnostics.py b/tools/read_insitu_diagnostics.py index c5425e204e..ceebcfa919 100755 --- a/tools/read_insitu_diagnostics.py +++ b/tools/read_insitu_diagnostics.py @@ -97,6 +97,42 @@ def energy_spread_eV(all_data, per_slice=False): else: return (constants.c**2 / constants.e) * gamma_spread(all_data["average"]) * all_data["mass"] +def temperature_in_ev_x(all_data, per_slice=True): + if per_slice: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_x(all_data)**2 * np.atleast_2d(all_data["mass"]).T + else: + return (constants.c**2 / constants.e) * normalized_momentum_std_x(all_data)**2 * np.atleast_2d(all_data["mass"]).T + else: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_x(all_data["average"])**2 * all_data["mass"] + else: + return (constants.c**2 / constants.e) * normalized_momentum_std_x(all_data["average"])**2 * all_data["mass"] + +def temperature_in_ev_y(all_data, per_slice=True): + if per_slice: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_y(all_data)**2 * np.atleast_2d(all_data["mass"]).T + else: + return (constants.c**2 / constants.e) * normalized_momentum_std_y(all_data)**2 * np.atleast_2d(all_data["mass"]).T + else: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_y(all_data["average"])**2 * all_data["mass"] + else: + return (constants.c**2 / constants.e) * normalized_momentum_std_y(all_data["average"])**2 * all_data["mass"] + +def temperature_in_ev_z(all_data, per_slice=True): + if per_slice: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_z(all_data)**2 * np.atleast_2d(all_data["mass"]).T + else: + return (constants.c**2 / constants.e) * normalized_momentum_std_z(all_data)**2 * np.atleast_2d(all_data["mass"]).T + else: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_z(all_data["average"])**2 * all_data["mass"] + else: + return (constants.c**2 / constants.e) * normalized_momentum_std_z(all_data["average"])**2 * all_data["mass"] + def position_std_x(all_data): return np.sqrt(np.maximum(all_data["[x^2]"] - all_data["[x]"]**2,0)) From d92e2bb371d4a59ba022896b2278e8d0066eddb4 Mon Sep 17 00:00:00 2001 From: Severin Date: Fri, 2 Jun 2023 22:47:00 +0200 Subject: [PATCH 05/11] remove typo --- src/particles/plasma/PlasmaParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index ef575d7a01..dbd04c910e 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -403,7 +403,7 @@ IonizationModule (const int lev, if(p_ion_mask[ip] != 0) { const long pid = amrex::Gpu::Atomic::Add( p_ip_elec, 1u ); - const long pidx = pid + old_size;beam + const long pidx = pid + old_size; // Copy ion data to new electron int_arrdata_elec[PlasmaIdx::id ][pidx] = pid_start + pid; From f0f0a235a6df248424e7f586712033f08136d821 Mon Sep 17 00:00:00 2001 From: Severin Date: Sat, 3 Jun 2023 11:35:12 +0200 Subject: [PATCH 06/11] add OpenMP support --- docs/source/run/parameters.rst | 4 +- .../plasma/PlasmaParticleContainer.cpp | 99 ++++++++++--------- 2 files changed, 52 insertions(+), 51 deletions(-) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index 686dab9ff6..846d0340b2 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -746,8 +746,8 @@ electron charge and mass. ``z_lo`` and ``z_hi`` are the lower and upper bounds of the z-axis of the simulation domain specified in the input file and can be used to generate a z/zeta-axis for plotting. ``normalized_density_factor`` is equal to ``dx * dy * dz`` in normalized units and 1 in -SI units. It can be used to convert ``sum(w)``, which specifies the beam density in normalized -units and beam weight an SI units, to the beam weight in both unit systems. +SI units. It can be used to convert ``sum(w)``, which specifies the particle density in normalized +units and particle weight in SI units, to the particle weight in both unit systems. The data is written to a file at ``/reduced_..txt``. The in-situ diagnostics file format consists of a header part in ASCII containing a JSON object. diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index dbd04c910e..ffd3dc8213 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -463,19 +463,18 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) // loading the data const auto ptd = pti.GetParticleTile().getParticleTileData(); - const int num_particles = pti.numParticles(); - // #ifdef AMREX_USE_OMP - // #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) - // #endif - // { - // amrex::Long const num_particles = pti.numParticles(); - // #ifdef AMREX_USE_OMP - // amrex::Long const idx_begin = (num_particles * omp_get_thread_num()) / omp_get_num_threads(); - // amrex::Long const idx_end = (num_particles * (omp_get_thread_num()+1)) / omp_get_num_threads(); - // #else - // amrex::Long constexpr idx_begin = 0; - // amrex::Long const idx_end = num_particles; - // #endif +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + { + amrex::Long const num_particles = pti.numParticles(); +#ifdef AMREX_USE_OMP + amrex::Long const idx_begin = (num_particles * omp_get_thread_num()) / omp_get_num_threads(); + amrex::Long const idx_end = (num_particles * (omp_get_thread_num()+1)) / omp_get_num_threads(); +#else + amrex::Long constexpr idx_begin = 0; + amrex::Long const idx_end = num_particles; +#endif // Tuple contains: // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 // sum(w), , , , , , , , , , , , , np @@ -490,10 +489,10 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) amrex::Real, int> reduce_data(reduce_op); using ReduceTuple = typename decltype(reduce_data)::Type; reduce_op.eval( - num_particles, reduce_data, - [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple + idx_end-idx_begin, reduce_data, + [=] AMREX_GPU_DEVICE (int idx) -> ReduceTuple { - const int ip = i; //indices[cell_start+i]; + const int ip = idx + idx_begin; const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); @@ -528,39 +527,41 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) 1}; }); - ReduceTuple a = reduce_data.value(); - const amrex::Real sum_w0 = amrex::get< 0>(a); - const amrex::Real sum_w_inv = sum_w0::epsilon() ? 0._rt : 1._rt/sum_w0; - - m_insitu_rdata[islice ] = sum_w0; - m_insitu_rdata[islice+ 1*m_nslices] = amrex::get< 1>(a)*sum_w_inv; - m_insitu_rdata[islice+ 2*m_nslices] = amrex::get< 2>(a)*sum_w_inv; - m_insitu_rdata[islice+ 3*m_nslices] = amrex::get< 3>(a)*sum_w_inv; - m_insitu_rdata[islice+ 4*m_nslices] = amrex::get< 4>(a)*sum_w_inv; - m_insitu_rdata[islice+ 5*m_nslices] = amrex::get< 5>(a)*sum_w_inv; - m_insitu_rdata[islice+ 6*m_nslices] = amrex::get< 6>(a)*sum_w_inv; - m_insitu_rdata[islice+ 7*m_nslices] = amrex::get< 7>(a)*sum_w_inv; - m_insitu_rdata[islice+ 8*m_nslices] = amrex::get< 8>(a)*sum_w_inv; - m_insitu_rdata[islice+ 9*m_nslices] = amrex::get< 9>(a)*sum_w_inv; - m_insitu_rdata[islice+10*m_nslices] = amrex::get<10>(a)*sum_w_inv; - m_insitu_rdata[islice+11*m_nslices] = amrex::get<11>(a)*sum_w_inv; - m_insitu_rdata[islice+12*m_nslices] = amrex::get<12>(a)*sum_w_inv; - m_insitu_idata[islice ] = amrex::get<13>(a); - - m_insitu_sum_rdata[ 0] += sum_w0; - m_insitu_sum_rdata[ 1] += amrex::get< 1>(a); - m_insitu_sum_rdata[ 2] += amrex::get< 2>(a); - m_insitu_sum_rdata[ 3] += amrex::get< 3>(a); - m_insitu_sum_rdata[ 4] += amrex::get< 4>(a); - m_insitu_sum_rdata[ 5] += amrex::get< 5>(a); - m_insitu_sum_rdata[ 6] += amrex::get< 6>(a); - m_insitu_sum_rdata[ 7] += amrex::get< 7>(a); - m_insitu_sum_rdata[ 8] += amrex::get< 8>(a); - m_insitu_sum_rdata[ 9] += amrex::get< 9>(a); - m_insitu_sum_rdata[10] += amrex::get<10>(a); - m_insitu_sum_rdata[11] += amrex::get<11>(a); - m_insitu_sum_rdata[12] += amrex::get<12>(a); - m_insitu_sum_idata[ 0] += amrex::get<13>(a); + ReduceTuple a = reduce_data.value(); + const amrex::Real sum_w0 = amrex::get< 0>(a); + const amrex::Real sum_w_inv = sum_w0::epsilon() ? 0._rt + : 1._rt/sum_w0; + + m_insitu_rdata[islice ] = sum_w0; + m_insitu_rdata[islice+ 1*m_nslices] = amrex::get< 1>(a)*sum_w_inv; + m_insitu_rdata[islice+ 2*m_nslices] = amrex::get< 2>(a)*sum_w_inv; + m_insitu_rdata[islice+ 3*m_nslices] = amrex::get< 3>(a)*sum_w_inv; + m_insitu_rdata[islice+ 4*m_nslices] = amrex::get< 4>(a)*sum_w_inv; + m_insitu_rdata[islice+ 5*m_nslices] = amrex::get< 5>(a)*sum_w_inv; + m_insitu_rdata[islice+ 6*m_nslices] = amrex::get< 6>(a)*sum_w_inv; + m_insitu_rdata[islice+ 7*m_nslices] = amrex::get< 7>(a)*sum_w_inv; + m_insitu_rdata[islice+ 8*m_nslices] = amrex::get< 8>(a)*sum_w_inv; + m_insitu_rdata[islice+ 9*m_nslices] = amrex::get< 9>(a)*sum_w_inv; + m_insitu_rdata[islice+10*m_nslices] = amrex::get<10>(a)*sum_w_inv; + m_insitu_rdata[islice+11*m_nslices] = amrex::get<11>(a)*sum_w_inv; + m_insitu_rdata[islice+12*m_nslices] = amrex::get<12>(a)*sum_w_inv; + m_insitu_idata[islice ] = amrex::get<13>(a); + + m_insitu_sum_rdata[ 0] += sum_w0; + m_insitu_sum_rdata[ 1] += amrex::get< 1>(a); + m_insitu_sum_rdata[ 2] += amrex::get< 2>(a); + m_insitu_sum_rdata[ 3] += amrex::get< 3>(a); + m_insitu_sum_rdata[ 4] += amrex::get< 4>(a); + m_insitu_sum_rdata[ 5] += amrex::get< 5>(a); + m_insitu_sum_rdata[ 6] += amrex::get< 6>(a); + m_insitu_sum_rdata[ 7] += amrex::get< 7>(a); + m_insitu_sum_rdata[ 8] += amrex::get< 8>(a); + m_insitu_sum_rdata[ 9] += amrex::get< 9>(a); + m_insitu_sum_rdata[10] += amrex::get<10>(a); + m_insitu_sum_rdata[11] += amrex::get<11>(a); + m_insitu_sum_rdata[12] += amrex::get<12>(a); + m_insitu_sum_idata[ 0] += amrex::get<13>(a); + } } } From 203cf0993f949a7f0d0c831460a9836d9513a319 Mon Sep 17 00:00:00 2001 From: Severin Date: Sat, 3 Jun 2023 16:05:06 +0200 Subject: [PATCH 07/11] remove omp support --- .../plasma/PlasmaParticleContainer.cpp | 88 ++++++++----------- 1 file changed, 37 insertions(+), 51 deletions(-) diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index ffd3dc8213..ea55bf15d0 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -463,18 +463,8 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) // loading the data const auto ptd = pti.GetParticleTile().getParticleTileData(); -#ifdef AMREX_USE_OMP -#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) -#endif - { - amrex::Long const num_particles = pti.numParticles(); -#ifdef AMREX_USE_OMP - amrex::Long const idx_begin = (num_particles * omp_get_thread_num()) / omp_get_num_threads(); - amrex::Long const idx_end = (num_particles * (omp_get_thread_num()+1)) / omp_get_num_threads(); -#else - amrex::Long constexpr idx_begin = 0; - amrex::Long const idx_end = num_particles; -#endif + amrex::Long const num_particles = pti.numParticles(); + // Tuple contains: // 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13 // sum(w), , , , , , , , , , , , , np @@ -489,11 +479,9 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) amrex::Real, int> reduce_data(reduce_op); using ReduceTuple = typename decltype(reduce_data)::Type; reduce_op.eval( - idx_end-idx_begin, reduce_data, - [=] AMREX_GPU_DEVICE (int idx) -> ReduceTuple + num_particles, reduce_data, + [=] AMREX_GPU_DEVICE (int ip) -> ReduceTuple { - const int ip = idx + idx_begin; - const amrex::Real xp = ptd.pos(0, ip); const amrex::Real yp = ptd.pos(1, ip); const amrex::Real uxp = ptd.rdata(PlasmaIdx::ux )[ip] * clight_inv; // proper velocity to u @@ -527,41 +515,39 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) 1}; }); - ReduceTuple a = reduce_data.value(); - const amrex::Real sum_w0 = amrex::get< 0>(a); - const amrex::Real sum_w_inv = sum_w0::epsilon() ? 0._rt - : 1._rt/sum_w0; - - m_insitu_rdata[islice ] = sum_w0; - m_insitu_rdata[islice+ 1*m_nslices] = amrex::get< 1>(a)*sum_w_inv; - m_insitu_rdata[islice+ 2*m_nslices] = amrex::get< 2>(a)*sum_w_inv; - m_insitu_rdata[islice+ 3*m_nslices] = amrex::get< 3>(a)*sum_w_inv; - m_insitu_rdata[islice+ 4*m_nslices] = amrex::get< 4>(a)*sum_w_inv; - m_insitu_rdata[islice+ 5*m_nslices] = amrex::get< 5>(a)*sum_w_inv; - m_insitu_rdata[islice+ 6*m_nslices] = amrex::get< 6>(a)*sum_w_inv; - m_insitu_rdata[islice+ 7*m_nslices] = amrex::get< 7>(a)*sum_w_inv; - m_insitu_rdata[islice+ 8*m_nslices] = amrex::get< 8>(a)*sum_w_inv; - m_insitu_rdata[islice+ 9*m_nslices] = amrex::get< 9>(a)*sum_w_inv; - m_insitu_rdata[islice+10*m_nslices] = amrex::get<10>(a)*sum_w_inv; - m_insitu_rdata[islice+11*m_nslices] = amrex::get<11>(a)*sum_w_inv; - m_insitu_rdata[islice+12*m_nslices] = amrex::get<12>(a)*sum_w_inv; - m_insitu_idata[islice ] = amrex::get<13>(a); - - m_insitu_sum_rdata[ 0] += sum_w0; - m_insitu_sum_rdata[ 1] += amrex::get< 1>(a); - m_insitu_sum_rdata[ 2] += amrex::get< 2>(a); - m_insitu_sum_rdata[ 3] += amrex::get< 3>(a); - m_insitu_sum_rdata[ 4] += amrex::get< 4>(a); - m_insitu_sum_rdata[ 5] += amrex::get< 5>(a); - m_insitu_sum_rdata[ 6] += amrex::get< 6>(a); - m_insitu_sum_rdata[ 7] += amrex::get< 7>(a); - m_insitu_sum_rdata[ 8] += amrex::get< 8>(a); - m_insitu_sum_rdata[ 9] += amrex::get< 9>(a); - m_insitu_sum_rdata[10] += amrex::get<10>(a); - m_insitu_sum_rdata[11] += amrex::get<11>(a); - m_insitu_sum_rdata[12] += amrex::get<12>(a); - m_insitu_sum_idata[ 0] += amrex::get<13>(a); - } + ReduceTuple a = reduce_data.value(); + const amrex::Real sum_w0 = amrex::get< 0>(a); + const amrex::Real sum_w_inv = sum_w0 <= 0._rt ? 0._rt : 1._rt/sum_w0; + + m_insitu_rdata[islice ] = sum_w0; + m_insitu_rdata[islice+ 1*m_nslices] = amrex::get< 1>(a)*sum_w_inv; + m_insitu_rdata[islice+ 2*m_nslices] = amrex::get< 2>(a)*sum_w_inv; + m_insitu_rdata[islice+ 3*m_nslices] = amrex::get< 3>(a)*sum_w_inv; + m_insitu_rdata[islice+ 4*m_nslices] = amrex::get< 4>(a)*sum_w_inv; + m_insitu_rdata[islice+ 5*m_nslices] = amrex::get< 5>(a)*sum_w_inv; + m_insitu_rdata[islice+ 6*m_nslices] = amrex::get< 6>(a)*sum_w_inv; + m_insitu_rdata[islice+ 7*m_nslices] = amrex::get< 7>(a)*sum_w_inv; + m_insitu_rdata[islice+ 8*m_nslices] = amrex::get< 8>(a)*sum_w_inv; + m_insitu_rdata[islice+ 9*m_nslices] = amrex::get< 9>(a)*sum_w_inv; + m_insitu_rdata[islice+10*m_nslices] = amrex::get<10>(a)*sum_w_inv; + m_insitu_rdata[islice+11*m_nslices] = amrex::get<11>(a)*sum_w_inv; + m_insitu_rdata[islice+12*m_nslices] = amrex::get<12>(a)*sum_w_inv; + m_insitu_idata[islice ] = amrex::get<13>(a); + + m_insitu_sum_rdata[ 0] += sum_w0; + m_insitu_sum_rdata[ 1] += amrex::get< 1>(a); + m_insitu_sum_rdata[ 2] += amrex::get< 2>(a); + m_insitu_sum_rdata[ 3] += amrex::get< 3>(a); + m_insitu_sum_rdata[ 4] += amrex::get< 4>(a); + m_insitu_sum_rdata[ 5] += amrex::get< 5>(a); + m_insitu_sum_rdata[ 6] += amrex::get< 6>(a); + m_insitu_sum_rdata[ 7] += amrex::get< 7>(a); + m_insitu_sum_rdata[ 8] += amrex::get< 8>(a); + m_insitu_sum_rdata[ 9] += amrex::get< 9>(a); + m_insitu_sum_rdata[10] += amrex::get<10>(a); + m_insitu_sum_rdata[11] += amrex::get<11>(a); + m_insitu_sum_rdata[12] += amrex::get<12>(a); + m_insitu_sum_idata[ 0] += amrex::get<13>(a); } } From bb873bea6dc6e02917e85d1239e39498cba75e36 Mon Sep 17 00:00:00 2001 From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com> Date: Tue, 6 Jun 2023 10:29:45 +0200 Subject: [PATCH 08/11] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Maxence Thévenet --- docs/source/run/parameters.rst | 23 +++++++++---------- .../plasma/PlasmaParticleContainer.H | 4 ++-- .../plasma/PlasmaParticleContainer.cpp | 4 ++-- 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst index 846d0340b2..de2831331a 100644 --- a/docs/source/run/parameters.rst +++ b/docs/source/run/parameters.rst @@ -722,10 +722,9 @@ Field diagnostics In-situ diagnostics ^^^^^^^^^^^^^^^^^^^ -Besides the standard diagnostics, fast in-situ diagnostics are available. They are most useful when beams with a large number of particles is used, because then the important moments can be calculated by GPU and significantly reduce the post-analysis of the simulation. -For particle beams, they can be used to calculate the main characterizing beam parameters (width, energy spread, emittance, etc.). Additionally, the plasma particle properties (e.g, the temperature) can be calculated. - -The in-situ quantities are calculated per slice and together with the weights per slice, all important quantities can be calculated. +Besides the standard diagnostics, fast in-situ diagnostics are available. They are most useful when beams with large numbers of particles are used, as the important moments can be calculated in-situ (during the simulation) to largely reduce the simulation's analysis. +In-situ diagnostics compute slice quantities (1 number per quantity per longitudinal cell). +For particle beams, they can be used to calculate the main characterizing beam parameters (width, energy spread, emittance, etc.), from which most common beam parameters (e.g. slice and projected emittance, etc.) can be computed. Additionally, the plasma particle properties (e.g, the temperature) can be calculated. For particle beams, the following quantities are calculated per slice and stored: ``sum(w), [x], [x^2], [y], [y^2], [z], [z^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], [x*ux], [y*uy], [z*uz], [ga], [ga^2], np``. @@ -740,11 +739,11 @@ Additionally, some metadata is also available: ``time, step, n_slices, charge, mass, z_lo, z_hi, normalized_density_factor``. ``time`` and ``step`` refers to the physical time of the simulation and step number of the current timestep. -``n_slices`` equal to the number of slices in the zeta direction. +``n_slices`` is the number of slices in the zeta direction. ``charge`` and ``mass`` relate to a single particle and are for example equal to the electron charge and mass. ``z_lo`` and ``z_hi`` are the lower and upper bounds of the z-axis of the simulation domain -specified in the input file and can be used to generate a z/zeta-axis for plotting. +specified in the input file and can be used to generate a z/zeta-axis for plotting (note that they corresponds to mesh nodes, while the data is cell-centered). ``normalized_density_factor`` is equal to ``dx * dy * dz`` in normalized units and 1 in SI units. It can be used to convert ``sum(w)``, which specifies the particle density in normalized units and particle weight in SI units, to the particle weight in both unit systems. @@ -752,19 +751,19 @@ units and particle weight in SI units, to the particle weight in both unit syste The data is written to a file at ``/reduced_..txt``. The in-situ diagnostics file format consists of a header part in ASCII containing a JSON object. When this is parsed into Python it can be converted to a NumPy structured datatype. -The rest of the file, following immediately after the closing }, is in binary format and -contains all of the in-situ diagnostic along with some meta data. This part can be read using the +The rest of the file, following immediately after the closing ``}``, is in binary format and +contains all of the in-situ diagnostics along with some metadata. This part can be read using the structured datatype of the first section. Use ``hipace/tools/read_insitu_diagnostics.py`` to read the files using this format. Functions to calculate the most useful properties are also provided in that file. -* `` or beams.insitu_period`` (`int`) optional (default ``-1``) - Period of the beam in-situ diagnostics. `-1` means no in-situ diagnostics. +* `` or beams.insitu_period`` (`int`) optional (default ``0``) + Period of the beam in-situ diagnostics. `0` means no beam in-situ diagnostics. * `` or beams.insitu_file_prefix`` (`string`) optional (default ``"diags/insitu"``) Path of the beam in-situ output. Must not be the same as `hipace.file_prefix`. -* `` or plasmas.insitu_period`` (`int`) optional (default ``-1``) - Period of the plasma in-situ diagnostics. `-1` means no in-situ diagnostics. +* `` or plasmas.insitu_period`` (`int`) optional (default ``0``) + Period of the plasma in-situ diagnostics. `0` means no plasma in-situ diagnostics. * `` or plasmas.insitu_file_prefix`` (`string`) optional (default ``"plasma_diags/insitu"``) Path of the plasma in-situ output. Must not be the same as `hipace.file_prefix`. diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index 3d2d3add6c..d4e1834b4f 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -138,7 +138,7 @@ public: /** Returns name of the plasma */ const std::string& GetName () const {return m_name;} - /** Compute reduced plasma diagnostics of current slice, store in member variable + /** Compute in-situ plasma diagnostics of current slice, store in member variable * \param[in] islice current slice, on which diags are computed. */ void InSituComputeDiags (int islice); @@ -215,7 +215,7 @@ private: /** Prefix/path for the output files */ std::string m_insitu_file_prefix = "diags/plasma_insitu"; /** How often the insitu plasma diagnostics should be computed and written */ - int m_insitu_period {-1}; + int m_insitu_period {0}; }; /** \brief Iterator over boxes in a particle container */ diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index ea55bf15d0..b6e12fe4ed 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -596,7 +596,7 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am {"[ux^2]" , &m_insitu_rdata[6*nslices], nslices}, {"[uy]" , &m_insitu_rdata[7*nslices], nslices}, {"[uy^2]" , &m_insitu_rdata[8*nslices], nslices}, - {"[uz]" , &m_insitu_rdata[9*nslices], nslices}, + {"[uz]" , &m_insitu_rdata[9*nslices], nslices}, {"[uz^2]" , &m_insitu_rdata[10*nslices], nslices}, {"[ga]" , &m_insitu_rdata[11*nslices], nslices}, {"[ga^2]" , &m_insitu_rdata[12*nslices], nslices}, @@ -611,7 +611,7 @@ PlasmaParticleContainer::InSituWriteToFile (int step, amrex::Real time, const am {"[ux^2]", &(m_insitu_sum_rdata[ 6] /= sum_w0)}, {"[uy]" , &(m_insitu_sum_rdata[ 7] /= sum_w0)}, {"[uy^2]", &(m_insitu_sum_rdata[ 8] /= sum_w0)}, - {"[uz]", &(m_insitu_sum_rdata[ 9] /= sum_w0)}, + {"[uz]" , &(m_insitu_sum_rdata[ 9] /= sum_w0)}, {"[uz^2]", &(m_insitu_sum_rdata[10] /= sum_w0)}, {"[ga]" , &(m_insitu_sum_rdata[11] /= sum_w0)}, {"[ga^2]", &(m_insitu_sum_rdata[12] /= sum_w0)} From 990993e72b8200cd77945aac778c9e7064980717 Mon Sep 17 00:00:00 2001 From: Severin Date: Tue, 6 Jun 2023 12:10:28 +0200 Subject: [PATCH 09/11] cleaning --- src/particles/beam/BeamParticleContainer.H | 10 +- src/particles/beam/BeamParticleContainer.cpp | 5 +- .../plasma/PlasmaParticleContainer.H | 4 +- .../plasma/PlasmaParticleContainer.cpp | 5 +- tools/read_insitu_diagnostics.py | 131 ++++++++++++------ 5 files changed, 99 insertions(+), 56 deletions(-) diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H index 14560e4112..8c7a408972 100644 --- a/src/particles/beam/BeamParticleContainer.H +++ b/src/particles/beam/BeamParticleContainer.H @@ -227,14 +227,13 @@ private: * 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, * sum(w), [x], [x^2], [y], [y^2], [z], [z^2], [ux], [ux^2], [uy], [uy^2], [uz], [uz^2], * - * 13, 14, 15, 16, 17, 18 - * [x*ux], [y*uy], [z*uz], [ga], [ga^2], np + * 13, 14, 15, 16, 17 + * [x*ux], [y*uy], [z*uz], [ga], [ga^2] * where [] means average over all particles within slice. * Per-slice emittance: sqrt( abs( ([x^2]-[x]^2) * ([ux^2]-[ux]^2) - ([x*ux]-[x][ux])^2 ) ). * Projected emittance: Same as above AFTER averaging all these quantities over slices. * Energy spread: sqrt([ga^2]-[ga]^2), and same as above. * Momentum spread: sqrt([uz^2]-[uz]^2), and same as above. - * Np: number of particles in this slice */ amrex::Vector m_insitu_rdata; /** Per-slice int beam properties: @@ -249,8 +248,9 @@ private: amrex::Vector m_insitu_sum_idata; /** Prefix/path for the output files */ std::string m_insitu_file_prefix = "diags/insitu"; - /** How often the insitu beam diagnostics should be computed and written */ - int m_insitu_period {-1}; + /** How often the insitu beam diagnostics should be computed and written + * Default is 0, meaning no output */ + int m_insitu_period {0}; }; #endif diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp index 71c472f81e..90b7d58a23 100644 --- a/src/particles/beam/BeamParticleContainer.cpp +++ b/src/particles/beam/BeamParticleContainer.cpp @@ -308,8 +308,7 @@ void BeamParticleContainer::TagByLevel (const int current_N_level, bool BeamParticleContainer::doInSitu (int step) { - if (m_insitu_period <= 0) return false; - return step % m_insitu_period == 0; + return (m_insitu_period > 0 && step % m_insitu_period == 0); } void @@ -322,7 +321,7 @@ BeamParticleContainer::InSituComputeDiags (int islice, int islice_local) AMREX_ALWAYS_ASSERT(m_insitu_rdata.size()>0 && m_insitu_idata.size()>0 && m_insitu_sum_rdata.size()>0 && m_insitu_sum_idata.size()>0); - PhysConst const phys_const = get_phys_const(); + const PhysConst phys_const = get_phys_const(); const amrex::Real clight_inv = 1.0_rt/phys_const.c; const amrex::Real clightsq_inv = 1.0_rt/(phys_const.c*phys_const.c); diff --git a/src/particles/plasma/PlasmaParticleContainer.H b/src/particles/plasma/PlasmaParticleContainer.H index d4e1834b4f..c8297b3df7 100644 --- a/src/particles/plasma/PlasmaParticleContainer.H +++ b/src/particles/plasma/PlasmaParticleContainer.H @@ -199,7 +199,6 @@ private: * Per-slice emittance: sqrt( abs( ([x^2]-[x]^2) * ([ux^2]-[ux]^2) - ([x*ux]-[x][ux])^2 ) ). * Projected emittance: Same as above AFTER averaging all these quantities over slices. * Energy spread: sqrt([ga^2]-[ga]^2), and same as above. - * Np: number of particles in this slice */ amrex::Vector m_insitu_rdata; /** Per-slice int plasma properties: @@ -214,7 +213,8 @@ private: amrex::Vector m_insitu_sum_idata; /** Prefix/path for the output files */ std::string m_insitu_file_prefix = "diags/plasma_insitu"; - /** How often the insitu plasma diagnostics should be computed and written */ + /** How often the insitu plasma diagnostics should be computed and written + * Default is 0, meaning no output */ int m_insitu_period {0}; }; diff --git a/src/particles/plasma/PlasmaParticleContainer.cpp b/src/particles/plasma/PlasmaParticleContainer.cpp index b6e12fe4ed..40b956db9d 100644 --- a/src/particles/plasma/PlasmaParticleContainer.cpp +++ b/src/particles/plasma/PlasmaParticleContainer.cpp @@ -439,8 +439,7 @@ IonizationModule (const int lev, bool PlasmaParticleContainer::doInSitu (int step) { - if (m_insitu_period <= 0) return false; - return step % m_insitu_period == 0; + return (m_insitu_period > 0 && step % m_insitu_period == 0); } void @@ -453,7 +452,7 @@ PlasmaParticleContainer::InSituComputeDiags (int islice) AMREX_ALWAYS_ASSERT(m_insitu_rdata.size()>0 && m_insitu_idata.size()>0 && m_insitu_sum_rdata.size()>0 && m_insitu_sum_idata.size()>0); - PhysConst const phys_const = get_phys_const(); + const PhysConst phys_const = get_phys_const(); const amrex::Real clight_inv = 1.0_rt/phys_const.c; const amrex::Real clightsq_inv = 1.0_rt/(phys_const.c*phys_const.c); diff --git a/tools/read_insitu_diagnostics.py b/tools/read_insitu_diagnostics.py index ceebcfa919..581839a39e 100755 --- a/tools/read_insitu_diagnostics.py +++ b/tools/read_insitu_diagnostics.py @@ -97,59 +97,104 @@ def energy_spread_eV(all_data, per_slice=False): else: return (constants.c**2 / constants.e) * gamma_spread(all_data["average"]) * all_data["mass"] -def temperature_in_ev_x(all_data, per_slice=True): - if per_slice: - if all_data["is_normalized_units"][0]: - return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_x(all_data)**2 * np.atleast_2d(all_data["mass"]).T - else: - return (constants.c**2 / constants.e) * normalized_momentum_std_x(all_data)**2 * np.atleast_2d(all_data["mass"]).T - else: - if all_data["is_normalized_units"][0]: - return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_x(all_data["average"])**2 * all_data["mass"] - else: - return (constants.c**2 / constants.e) * normalized_momentum_std_x(all_data["average"])**2 * all_data["mass"] +def temperature_in_ev(all_data, per_slice=True, direction='all'): + """ + calculated temperature in eV for in-situ diagnostics -def temperature_in_ev_y(all_data, per_slice=True): - if per_slice: - if all_data["is_normalized_units"][0]: - return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_y(all_data)**2 * np.atleast_2d(all_data["mass"]).T - else: - return (constants.c**2 / constants.e) * normalized_momentum_std_y(all_data)**2 * np.atleast_2d(all_data["mass"]).T - else: - if all_data["is_normalized_units"][0]: - return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_y(all_data["average"])**2 * all_data["mass"] - else: - return (constants.c**2 / constants.e) * normalized_momentum_std_y(all_data["average"])**2 * all_data["mass"] + Parameters + ---------- -def temperature_in_ev_z(all_data, per_slice=True): - if per_slice: - if all_data["is_normalized_units"][0]: - return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_z(all_data)**2 * np.atleast_2d(all_data["mass"]).T + all_data: NumPy structured array over timesteps. + To be obtained via read_file() (see above) + per_slice: boolean, whether the temperature should be returned per slice or averaged over all slices + direction: direction along which the temperature should be calculated. Available options: + 'x', 'y', 'z', 'all' + + Returns + ------- + + NumPy structured array over timesteps with the tempature in eV. + """ + if direction=='all': + return 1/3.*(temperature_in_ev(all_data, per_slice=per_slice, direction='x') + + temperature_in_ev(all_data, per_slice=per_slice, direction='y') + + temperature_in_ev(all_data, per_slice=per_slice, direction='z')) + elif (direction=='x' or direction=='y' or direction=='z': + if per_slice: + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std(all_data, direction=direction)**2 * np.atleast_2d(all_data["mass"]).T + else: + return (constants.c**2 / constants.e) * normalized_momentum_std(all_data, direction=direction)**2 * np.atleast_2d(all_data["mass"]).T else: - return (constants.c**2 / constants.e) * normalized_momentum_std_z(all_data)**2 * np.atleast_2d(all_data["mass"]).T + if all_data["is_normalized_units"][0]: + return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std(all_data["average"], direction=direction)**2 * all_data["mass"] + else: + return (constants.c**2 / constants.e) * normalized_momentum_std(all_data["average"], direction=direction)**2 * all_data["mass"] else: - if all_data["is_normalized_units"][0]: - return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std_z(all_data["average"])**2 * all_data["mass"] - else: - return (constants.c**2 / constants.e) * normalized_momentum_std_z(all_data["average"])**2 * all_data["mass"] + print("Error, unknown direction, use 'x', 'y', 'z', or 'all'") + return + +def position_std(all_data, direction='x'): + """ + calculated standard deviation of the position in x,y, or z for in-situ diagnostics + + Parameters + ---------- + + all_data: NumPy structured array over timesteps. + To be obtained via read_file() (see above) + direction: direction along which the temperature should be calculated. Available options: + 'x', 'y', 'z' + + Returns + ------- -def position_std_x(all_data): - return np.sqrt(np.maximum(all_data["[x^2]"] - all_data["[x]"]**2,0)) + NumPy structured array over timesteps. + """ + if direction=='x': + var = "[x]" + var_sq = "[x^2]" + elif direction=='y': + var = "[y]" + var_sq = "[y^2]" + elif direction=='z': + var = "[z]" + var_sq = "[z^2]" + else: + print("Error, unknown direction, use 'x', 'y', or 'z'") + return + return np.sqrt(np.maximum(all_data[var_sq] - all_data[var]**2,0)) -def position_std_y(all_data): - return np.sqrt(np.maximum(all_data["[y^2]"] - all_data["[y]"]**2,0)) +def normalized_momentum_std(all_data, direction='x'): + """ + calculated standard deviation of the momentum in x,y, or z for in-situ diagnostics -def position_std_z(all_data): - return np.sqrt(np.maximum(all_data["[z^2]"] - all_data["[z]"]**2,0)) + Parameters + ---------- -def normalized_momentum_std_x(all_data): - return np.sqrt(np.maximum(all_data["[ux^2]"] - all_data["[ux]"]**2,0)) + all_data: NumPy structured array over timesteps. + To be obtained via read_file() (see above) + direction: direction along which the temperature should be calculated. Available options: + 'x', 'y', 'z' -def normalized_momentum_std_y(all_data): - return np.sqrt(np.maximum(all_data["[uy^2]"] - all_data["[uy]"]**2,0)) + Returns + ------- -def normalized_momentum_std_z(all_data): - return np.sqrt(np.maximum(all_data["[uz^2]"] - all_data["[uz]"]**2,0)) + NumPy structured array over timesteps. + """ + if direction=='x': + var = "[ux]" + var_sq = "[ux^2]" + elif direction=='y': + var = "[uy]" + var_sq = "[uy^2]" + elif direction=='z': + var = "[uz]" + var_sq = "[uz^2]" + else: + print("Error, unknown direction, use 'x', 'y', or 'z'") + return + return np.sqrt(np.maximum(all_data[var_sq] - all_data[var]**2,0)) def per_slice_charge(all_data): return np.atleast_2d(all_data["charge"]).T * all_data["sum(w)"] * np.atleast_2d(all_data["normalized_density_factor"]).T From 56d8b4dc0696627637a611df634e11ec65ecf15b Mon Sep 17 00:00:00 2001 From: Severin Diederichs <65728274+SeverinDiederichs@users.noreply.github.com> Date: Tue, 6 Jun 2023 12:37:07 +0200 Subject: [PATCH 10/11] Update tools/read_insitu_diagnostics.py Typo --- tools/read_insitu_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/read_insitu_diagnostics.py b/tools/read_insitu_diagnostics.py index 581839a39e..0d4371977f 100755 --- a/tools/read_insitu_diagnostics.py +++ b/tools/read_insitu_diagnostics.py @@ -119,7 +119,7 @@ def temperature_in_ev(all_data, per_slice=True, direction='all'): return 1/3.*(temperature_in_ev(all_data, per_slice=per_slice, direction='x') + temperature_in_ev(all_data, per_slice=per_slice, direction='y') + temperature_in_ev(all_data, per_slice=per_slice, direction='z')) - elif (direction=='x' or direction=='y' or direction=='z': + elif (direction=='x' or direction=='y' or direction=='z'): if per_slice: if all_data["is_normalized_units"][0]: return (constants.m_e * constants.c**2 / constants.e) * normalized_momentum_std(all_data, direction=direction)**2 * np.atleast_2d(all_data["mass"]).T From a870e949851777c312be35d7682db152fd308576 Mon Sep 17 00:00:00 2001 From: Severin Date: Wed, 7 Jun 2023 16:07:41 +0200 Subject: [PATCH 11/11] make eV consistent --- tools/read_insitu_diagnostics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tools/read_insitu_diagnostics.py b/tools/read_insitu_diagnostics.py index 0d4371977f..80f3418a65 100755 --- a/tools/read_insitu_diagnostics.py +++ b/tools/read_insitu_diagnostics.py @@ -97,7 +97,7 @@ def energy_spread_eV(all_data, per_slice=False): else: return (constants.c**2 / constants.e) * gamma_spread(all_data["average"]) * all_data["mass"] -def temperature_in_ev(all_data, per_slice=True, direction='all'): +def temperature_in_eV(all_data, per_slice=True, direction='all'): """ calculated temperature in eV for in-situ diagnostics @@ -116,9 +116,9 @@ def temperature_in_ev(all_data, per_slice=True, direction='all'): NumPy structured array over timesteps with the tempature in eV. """ if direction=='all': - return 1/3.*(temperature_in_ev(all_data, per_slice=per_slice, direction='x') + - temperature_in_ev(all_data, per_slice=per_slice, direction='y') + - temperature_in_ev(all_data, per_slice=per_slice, direction='z')) + return 1/3.*(temperature_in_eV(all_data, per_slice=per_slice, direction='x') + + temperature_in_eV(all_data, per_slice=per_slice, direction='y') + + temperature_in_eV(all_data, per_slice=per_slice, direction='z')) elif (direction=='x' or direction=='y' or direction=='z'): if per_slice: if all_data["is_normalized_units"][0]: