aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt1
-rw-r--r--Source/Diagnostics/ComputeDiagFunctors/Make.package1
-rw-r--r--Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H50
-rw-r--r--Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp130
-rw-r--r--Source/Diagnostics/Diagnostics.H15
-rw-r--r--Source/Diagnostics/Diagnostics.cpp85
-rw-r--r--Source/Diagnostics/FullDiagnostics.cpp18
7 files changed, 285 insertions, 15 deletions
diff --git a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt
index 5e085db4d..b2d0c096d 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt
+++ b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt
@@ -8,4 +8,5 @@ target_sources(WarpX
PartPerGridFunctor.cpp
BackTransformFunctor.cpp
BackTransformParticleFunctor.cpp
+ ParticleReductionFunctor.cpp
)
diff --git a/Source/Diagnostics/ComputeDiagFunctors/Make.package b/Source/Diagnostics/ComputeDiagFunctors/Make.package
index 7a07273e2..611ccb60d 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/Make.package
+++ b/Source/Diagnostics/ComputeDiagFunctors/Make.package
@@ -6,5 +6,6 @@ CEXE_sources += DivEFunctor.cpp
CEXE_sources += RhoFunctor.cpp
CEXE_sources += BackTransformFunctor.cpp
CEXE_sources += BackTransformParticleFunctor.cpp
+CEXE_sources += ParticleReductionFunctor.cpp
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics/ComputeDiagFunctors
diff --git a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H
new file mode 100644
index 000000000..81172a423
--- /dev/null
+++ b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.H
@@ -0,0 +1,50 @@
+#ifndef WARPX_PARTICLEREDUCTIONFUNCTOR_H_
+#define WARPX_PARTICLEREDUCTIONFUNCTOR_H_
+
+#include "ComputeDiagFunctor.H"
+
+#include <AMReX_Parser.H>
+
+#include <AMReX_BaseFwd.H>
+
+#include <memory>
+#include <string>
+
+/**
+ * \brief Functor to calculate per-cell averages of particle properties.
+ */
+class
+ParticleReductionFunctor final : public ComputeDiagFunctor
+{
+public:
+ /** Constructor.
+ * \param[in] mf_src source multifab. Must be nullptr as no source MF is needed
+ * to compute the number of particles per cell, banane.
+ * \param[in] lev level of multifab.
+ * \param[in] crse_ratio for interpolating field values from simulation MultiFabs
+ to the output diagnostic MultiFab mf_dst.
+ * \param[in] fn_str parser string that describes the function to apply to particles
+ * and then average over each cell
+ * \param[in] ispec index of the species over which to calculate the average
+ * \param[in] ncomp Number of component of mf_src to cell-center in dst multifab.
+ */
+ ParticleReductionFunctor(const amrex::MultiFab * const mf_src, const int lev,
+ const amrex::IntVect crse_ratio, const std::string fn_str, const int ispec, const int ncomp=1);
+
+ /** \brief Compute the average of the function m_map_fn over each grid cell.
+ *
+ * \param[out] mf_dst output MultiFab where the result is written
+ * \param[in] dcomp first component of mf_dst in which cell-centered
+ * data is stored
+ */
+ virtual void operator()(amrex::MultiFab& mf_dst, const int dcomp, const int /*i_buffer=0*/) const override;
+private:
+ int const m_lev; /**< level on which mf_src is defined */
+ int const m_ispec; /**< index of species to average over */
+ /** Parser function to be averaged by the functor. Arguments: x, y, z, ux, uy, uz */
+ std::unique_ptr<amrex::Parser> m_map_fn_parser;
+ /** Compiled #m_map_fn_parser */
+ amrex::ParserExecutor<6> m_map_fn;
+};
+
+#endif // WARPX_PARTICLEREDUCTIONFUNCTOR_H_
diff --git a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp
new file mode 100644
index 000000000..c5f9028f5
--- /dev/null
+++ b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp
@@ -0,0 +1,130 @@
+
+#include "ParticleReductionFunctor.H"
+
+#include "Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H"
+#include "Particles/MultiParticleContainer.H"
+#include "Particles/WarpXParticleContainer.H"
+#include "Utils/CoarsenIO.H"
+#include "WarpX.H"
+
+#include <AMReX_Array.H>
+#include <AMReX_BLassert.H>
+#include <AMReX_IntVect.H>
+#include <AMReX_MultiFab.H>
+#include <AMReX_REAL.H>
+
+using namespace amrex::literals;
+
+ParticleReductionFunctor::ParticleReductionFunctor (const amrex::MultiFab* mf_src, const int lev,
+ const amrex::IntVect crse_ratio, const std::string fn_str, const int ispec, const int ncomp)
+ : ComputeDiagFunctor(ncomp, crse_ratio), m_lev(lev), m_ispec(ispec)
+{
+ // mf_src will not be used, let's make sure it's null.
+ AMREX_ALWAYS_ASSERT(mf_src == nullptr);
+ // Write only in one output component.
+ AMREX_ALWAYS_ASSERT(ncomp == 1);
+
+ // Allocate and compile a parser based on the input string fn_str
+ m_map_fn_parser = std::make_unique<amrex::Parser>(makeParser(
+ fn_str, {"x", "y", "z", "ux", "uy", "uz"}));
+ m_map_fn = m_map_fn_parser->compile<6>();
+}
+
+void
+ParticleReductionFunctor::operator() (amrex::MultiFab& mf_dst, const int dcomp, const int /*i_buffer*/) const
+{
+ auto& warpx = WarpX::GetInstance();
+ // Guard cell is set to 1 for generality. However, for a cell-centered
+ // output Multifab, mf_dst, the guard-cell data is not needed especially considering
+ // the operations performend in the CoarsenAndInterpolate function.
+ constexpr int ng = 1;
+ // Temporary cell-centered, multi-component MultiFab for storing particles per cell.
+ amrex::MultiFab red_mf(warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), 1, ng);
+ amrex::MultiFab ppc_mf(warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), 1, ng);
+ // Set value to 0, and increment the value in each cell with ppc.
+ red_mf.setVal(0._rt);
+ ppc_mf.setVal(0._rt);
+ auto& pc = warpx.GetPartContainer().GetParticleContainer(m_ispec);
+ // Copy over compiled parser function so that it can be captured by value in the lambda
+ auto map_fn = m_map_fn;
+ ParticleToMesh(pc, red_mf, m_lev,
+ [=] AMREX_GPU_DEVICE (const WarpXParticleContainer::SuperParticleType& p,
+ amrex::Array4<amrex::Real> const& out_array,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi)
+ {
+ // Get position in WarpX convention to use in parser. Will be different from
+ // p.pos() for 1D and 2D simulations.
+ amrex::ParticleReal xw = 0._rt, yw = 0._rt, zw = 0._rt;
+ get_particle_position(p, xw, yw, zw);
+
+ // Get position in AMReX convention to calculate corresponding index.
+ // Ideally this will be replaced with the AMReX NGP interpolator
+ // Always do x direction. No RZ case because it's not implemented, and code
+ // will have aborted
+ int ii = 0, jj = 0, kk = 0;
+ amrex::ParticleReal x = p.pos(0);
+ amrex::Real lx = (x - plo[0]) * dxi[0];
+ ii = static_cast<int>(amrex::Math::floor(lx));
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_3D)
+ amrex::ParticleReal y = p.pos(1);
+ amrex::Real ly = (y - plo[1]) * dxi[1];
+ jj = static_cast<int>(amrex::Math::floor(ly));
+#endif
+#if defined(WARPX_DIM_3D)
+ amrex::ParticleReal z = p.pos(2);
+ amrex::Real lz = (z - plo[2]) * dxi[2];
+ kk = static_cast<int>(amrex::Math::floor(lz));
+#endif
+
+ // Fix dimensions since parser assumes u = gamma * v / c
+ amrex::ParticleReal ux = p.rdata(PIdx::ux) / PhysConst::c;
+ amrex::ParticleReal uy = p.rdata(PIdx::uy) / PhysConst::c;
+ amrex::ParticleReal uz = p.rdata(PIdx::uz) / PhysConst::c;
+ amrex::Real value = map_fn(xw, yw, zw, ux, uy, uz);
+ amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), p.rdata(PIdx::w) * value);
+ });
+ // Add the weight for each particle -- total number of particles of this species
+ ParticleToMesh(pc, ppc_mf, m_lev,
+ [=] AMREX_GPU_DEVICE (const WarpXParticleContainer::SuperParticleType& p,
+ amrex::Array4<amrex::Real> const& out_array,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& plo,
+ amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi)
+ {
+
+ // Get position in AMReX convention to calculate corresponding index.
+ // Ideally this will be replaced with the AMReX NGP interpolator
+ // Always do x direction. No RZ case because it's not implemented, and code
+ // will have aborted
+ int ii = 0, jj = 0, kk = 0;
+ amrex::ParticleReal x = p.pos(0);
+ amrex::Real lx = (x - plo[0]) * dxi[0];
+ ii = static_cast<int>(amrex::Math::floor(lx));
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_3D)
+ amrex::ParticleReal y = p.pos(1);
+ amrex::Real ly = (y - plo[1]) * dxi[1];
+ jj = static_cast<int>(amrex::Math::floor(ly));
+#endif
+#if defined(WARPX_DIM_3D)
+ amrex::ParticleReal z = p.pos(2);
+ amrex::Real lz = (z - plo[2]) * dxi[2];
+ kk = static_cast<int>(amrex::Math::floor(lz));
+#endif
+ amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), p.rdata(PIdx::w));
+ });
+ // Divide value by number of particles for average. Set average to zero if there are no particles
+ for (amrex::MFIter mfi(red_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
+ {
+ const amrex::Box& box = mfi.tilebox();
+ amrex::Array4<amrex::Real> const& a_red = red_mf.array(mfi);
+ amrex::Array4<amrex::Real const> const& a_ppc = ppc_mf.const_array(mfi);
+ amrex::ParallelFor(box,
+ [=] AMREX_GPU_DEVICE (int i, int j, int k) {
+ if (a_ppc(i,j,k,0) == 0) a_red(i,j,k,0) = 0;
+ else a_red(i,j,k,0) = a_red(i,j,k,0) / a_ppc(i,j,k,0);
+ });
+ }
+
+ // Coarsen and interpolate from ppc_mf to the output diagnostic MultiFab, mf_dst.
+ CoarsenIO::Coarsen(mf_dst, red_mf, dcomp, 0, nComp(), 0, m_crse_ratio);
+}
diff --git a/Source/Diagnostics/Diagnostics.H b/Source/Diagnostics/Diagnostics.H
index a88e6221c..4b90cdb22 100644
--- a/Source/Diagnostics/Diagnostics.H
+++ b/Source/Diagnostics/Diagnostics.H
@@ -134,9 +134,20 @@ protected:
/** Index of diagnostics in MultiDiagnostics::alldiags */
int m_diag_index;
/** Names of each component requested by the user.
- * in cylindrical geometry, this list is appended with
- * automatically-constructed names for all modes of all fields */
+ * The list is appended with the average particle fields, if used.
+ * In cylindrical geometry, this list is also appended with
+ * automatically-constructed names for all modes of all fields.*/
amrex::Vector< std::string > m_varnames;
+ /** Names of plotfile fields requested by the user */
+ amrex::Vector< std::string > m_varnames_fields;
+ /** Names of averaged particle properties to output */
+ amrex::Vector< std::string > m_pfield_varnames;
+ /** Names of species for which to output average particle fields */
+ std::vector< std::string > m_pfield_species;
+ /** Species indices corresponding to elements of m_pfield_varnames */
+ std::vector< int > m_pfield_species_index;
+ /** Map connecting particle field variable names to the corresponding parser strings */
+ std::map< std::string, std::string > m_pfield_strings = {};
/** If true, a dump is performed at the last timestep regardless of the required dump timesteps */
bool m_dump_last_timestep = true;
/** format for output files, "plotfile" or "openpmd" or "sensei" or "ascent"
diff --git a/Source/Diagnostics/Diagnostics.cpp b/Source/Diagnostics/Diagnostics.cpp
index fb41bebf1..b36011c1c 100644
--- a/Source/Diagnostics/Diagnostics.cpp
+++ b/Source/Diagnostics/Diagnostics.cpp
@@ -56,27 +56,27 @@ Diagnostics::BaseReadParameters ()
pp_diag_name.query("dump_last_timestep", m_dump_last_timestep);
// Query list of grid fields to write to output
- bool varnames_specified = pp_diag_name.queryarr("fields_to_plot", m_varnames);
+ bool varnames_specified = pp_diag_name.queryarr("fields_to_plot", m_varnames_fields);
if (!varnames_specified){
- m_varnames = {"Ex", "Ey", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz"};
+ m_varnames_fields = {"Ex", "Ey", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz"};
}
// Sanity check if user requests to plot phi
- if (WarpXUtilStr::is_in(m_varnames, "phi")){
+ if (WarpXUtilStr::is_in(m_varnames_fields, "phi")){
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
warpx.do_electrostatic==ElectrostaticSolverAlgo::LabFrame,
"plot phi only works if do_electrostatic = labframe");
}
// Sanity check if user requests to plot F
- if (WarpXUtilStr::is_in(m_varnames, "F")){
+ if (WarpXUtilStr::is_in(m_varnames_fields, "F")){
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
warpx.do_dive_cleaning,
"plot F only works if warpx.do_dive_cleaning = 1");
}
// G can be written to file only if WarpX::do_divb_cleaning = 1
- if (WarpXUtilStr::is_in(m_varnames, "G"))
+ if (WarpXUtilStr::is_in(m_varnames_fields, "G"))
{
WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
warpx.do_divb_cleaning, "G can be written to file only if warpx.do_divb_cleaning = 1");
@@ -85,9 +85,74 @@ Diagnostics::BaseReadParameters ()
// If user requests to plot proc_number for a serial run,
// delete proc_number from fields_to_plot
if (amrex::ParallelDescriptor::NProcs() == 1){
- m_varnames.erase(
- std::remove(m_varnames.begin(), m_varnames.end(), "proc_number"),
- m_varnames.end());
+ m_varnames_fields.erase(
+ std::remove(m_varnames_fields.begin(), m_varnames_fields.end(), "proc_number"),
+ m_varnames_fields.end());
+ }
+
+ // Get names of particle quantities to average at each grid point
+ const bool pfield_varnames_specified = pp_diag_name.queryarr("particle_fields_to_plot", m_pfield_varnames);
+ if (!pfield_varnames_specified){
+ m_pfield_varnames = {};
+ }
+#ifdef WARPX_DIM_RZ
+ if (m_pfield_varnames.size() != 0) {
+ amrex::Abort("Input error: cannot use particle_fields_to_plot not implemented for RZ");
+ }
+#endif
+
+
+ // Get parser strings for particle fields and generate map of parsers
+ std::string parser_str;
+ amrex::ParmParse pp_diag_pfield(m_diag_name + ".particle_fields");
+ for (const auto& var : m_pfield_varnames) {
+ Store_parserString(pp_diag_pfield, (var + "(x,y,z,ux,uy,uz)").c_str(), parser_str);
+ if (parser_str != "") {
+ m_pfield_strings.insert({var, parser_str});
+ }
+ else {
+ amrex::Abort("Input error: cannot find parser string for " + var + "." +
+ m_diag_name + ".particle_fields." + var + " in file");
+ }
+ }
+
+ // Names of all species in the simulation
+ m_all_species_names = warpx.GetPartContainer().GetSpeciesNames();
+
+ // Get names of species to average at each grid point
+ const bool pfield_species_specified = pp_diag_name.queryarr("particle_fields_species", m_pfield_species);
+ if (!pfield_species_specified){
+ m_pfield_species = m_all_species_names;
+ }
+
+ // Check that species names specified in m_pfield_species are valid
+ bool p_species_name_is_wrong;
+ // Loop over all species specified above
+ for (const auto& species : m_pfield_species) {
+ // Boolean used to check if species name was misspelled
+ p_species_name_is_wrong = true;
+ // Loop over all species
+ for (int i = 0, n = int(m_all_species_names.size()); i < n; i++) {
+ if (species == m_all_species_names[i]) {
+ // Store species index: will be used in ParticleReductionFunctor to calculate
+ // averages for this species
+ m_pfield_species_index.push_back(i);
+ p_species_name_is_wrong = false;
+ }
+ }
+ // If species name was misspelled, abort with error message
+ if (p_species_name_is_wrong) {
+ amrex::Abort("Input error: string " + species + " in " + m_diag_name +
+ ".particle_fields_species does not match any species");
+ }
+ }
+
+ m_varnames = m_varnames_fields;
+ // Generate names of averaged particle fields and append to m_varnames
+ for (int ivar=0; ivar<m_pfield_varnames.size(); ivar++) {
+ for (int ispec=0; ispec < int(m_pfield_species.size()); ispec++) {
+ m_varnames.push_back(m_pfield_varnames[ivar] + '_' + m_pfield_species[ispec]);
+ }
}
// Read user-defined physical extents for the output and store in m_lo and m_hi.
@@ -139,8 +204,6 @@ Diagnostics::BaseReadParameters ()
// Names of species to write to output
bool species_specified = pp_diag_name.queryarr("species", m_output_species_names);
- // Names of all species in the simulation
- m_all_species_names = warpx.GetPartContainer().GetSpeciesNames();
// Auxiliary variables
std::string species;
@@ -175,6 +238,8 @@ Diagnostics::BaseReadParameters ()
bool checkpoint_compatibility = false;
if (m_format == "checkpoint"){
if ( varnames_specified == false &&
+ pfield_varnames_specified == false &&
+ pfield_species_specified == false &&
lo_specified == false &&
hi_specified == false &&
cr_specified == false &&
diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp
index fb7a055f6..79791ae8e 100644
--- a/Source/Diagnostics/FullDiagnostics.cpp
+++ b/Source/Diagnostics/FullDiagnostics.cpp
@@ -5,6 +5,7 @@
#include "ComputeDiagFunctors/DivEFunctor.H"
#include "ComputeDiagFunctors/PartPerCellFunctor.H"
#include "ComputeDiagFunctors/PartPerGridFunctor.H"
+#include "ComputeDiagFunctors/ParticleReductionFunctor.H"
#include "ComputeDiagFunctors/RhoFunctor.H"
#include "Diagnostics/Diagnostics.H"
#include "Diagnostics/ParticleDiag/ParticleDiag.H"
@@ -408,10 +409,13 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
// Species index to loop over species that dump rho per species
int i = 0;
- m_all_field_functors[lev].resize( m_varnames.size() );
+ const auto nvar = static_cast<int>(m_varnames_fields.size());
+ const auto nspec = static_cast<int>(m_pfield_species.size());
+ const auto ntot = static_cast<int>(nvar + m_pfield_varnames.size() * nspec);
+
+ m_all_field_functors[lev].resize(ntot);
// Fill vector of functors for all components except individual cylindrical modes.
- const auto n = static_cast<int>(m_all_field_functors[lev].size());
- for (int comp=0; comp<n; comp++){
+ for (int comp=0; comp<nvar; comp++){
if ( m_varnames[comp] == "Ex" ){
m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Efield_aux(lev, 0), lev, m_crse_ratio);
} else if ( m_varnames[comp] == "Ey" ){
@@ -456,6 +460,14 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
amrex::Abort("Error: " + m_varnames[comp] + " is not a known field output type");
}
}
+ // Add functors for average particle data for each species
+ for (int pcomp=0; pcomp<int(m_pfield_varnames.size()); pcomp++) {
+ std::string varname = m_pfield_varnames[pcomp];
+ for (int ispec=0; ispec<int(m_pfield_species.size()); ispec++) {
+ m_all_field_functors[lev][nvar + pcomp * nspec + ispec] = std::make_unique<ParticleReductionFunctor>(nullptr,
+ lev, m_crse_ratio, m_pfield_strings[varname], m_pfield_species_index[ispec]);
+ }
+ }
AddRZModesToDiags( lev );
}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-18 05:17:15 +0000