diff options
author | 2022-01-14 09:30:14 -0800 | |
---|---|---|
committer | 2022-01-14 17:30:14 +0000 | |
commit | 2efde6cd40e0a155f188c4f53d96a11aa9ac13ac (patch) | |
tree | 72ae24c05533b366ec29ad2c387a4e787167dd11 /Source | |
parent | 400f536bed944f2e430f7c34dbb3fc00ac10332b (diff) | |
download | WarpX-2efde6cd40e0a155f188c4f53d96a11aa9ac13ac.tar.gz WarpX-2efde6cd40e0a155f188c4f53d96a11aa9ac13ac.tar.zst WarpX-2efde6cd40e0a155f188c4f53d96a11aa9ac13ac.zip |
Field probe line detector (#2513)
* FieldProbe using Particle
Update FieldProbe.cpp
Update FieldProbeParticleContainer.H
Updates FieldProbe and FieldProbeParticleContainer
* Make <diag>.integrate optional
The param parser query keeps te default value if no entry is found.
* Fixed number particle needed for AddNParticles
* Removing unnecessary type definition
* Added Doxygen-style comments to FieldProbe.cpp
Corrected Poynting calculation by implementing vacuum permeability
* Added Doxygen comments
* Implement virtual function ReducedDiags::AllocData() + comments
* InitData implemented
* Fixed Doxygen commenting.
* Now uses WarpX physics constant for vaccuum permeability
* forgotton comments to MultiReducedDiags
* Update Source/Diagnostics/ReducedDiags/FieldProbe.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update FieldProbe.H
* Update FieldProbe.cpp
* Update Source/Diagnostics/ReducedDiags/ReducedDiags.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/MultiReducedDiags.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/MultiReducedDiags.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update FieldProbeParticleContainer.H
* Update Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update FieldProbeParticleContainer.cpp
* Update FieldProbe.cpp
* Update FieldProbe.H
* Update Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.H
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update FieldProbeParticleContainer.cpp
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/ReducedDiags.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Changed enumerated class to struct w/ enumeration. Can remove "static_cast<int>"
* FieldProbeParticleContainer::iterator implemented
* Cleaned up output += operator, fixed output comments
* style fix
* Replaces Tabs with 4 spaces
* Defined modes and interp order to avoid GPU compilation errors
* 1 more tab fix
* EoL white spaces
* fixed a typoX
* Explicitly capturing "this" in parallel for to combat error saying "error #3223-D: Implicit capture of 'this' in extended lambda expression"
* removed unncessacesy double define
* moved output out of ParallelFor. temp variable for integrate
* Parse integrate, integrate all time steps, output setup for integrate and regular
* Fixed integrate bug.
* ammend header. integreate variable name change.
* Integrate values in input file
* updates to timing for integrate
* Update Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* whitespace
* Update Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Functionality to create 2D line of particles. Input included. No output yet
* ammend compiler errors
* Apply suggestions from code review - Style
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update reduce_diag_names
* 2D array setup- not complete
* field_probe_integrate change
* review amends
* Apply suggestions from code review - Style
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Vectors + AddNParticles
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update FieldProbe.cpp
* Update FieldProbe.H
* bug fix and inputs
* reintroduce raw_fields functionality
* docs update and correction
* whitespaces
* Fix GPU Compile (raw_fields)
* changed f_probe to m_probe apropriately
* Typos
Co-authored-by: David Grote <dpgrote@lbl.gov>
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Better name for ParticleVal
* used map for observables and units
* Apply suggestions from code review
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Simplified output. Fixed double integration error
* Update FieldProbe.H
Removed unneeded variable
* comments and fixed rawFields
* white spaces
* Apply suggestions from code review
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Update FieldProbe.H
* Guard on write
* Update Source/Diagnostics/ReducedDiags/FieldProbe.cpp
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Fix Syntax Error in Write
* Fix Init: Only 1 Particle (MPI)
Only one MPI rank adds a particle, which we done distribute
into the right rank.
* Fix MPI Deadlock: No Early Return
We just want to skip the write to `m_data`, not the rest
of the logic.
* Vector storage, Add N particle, debugging
* Fix Probe in Domain Logic
General global check, not only on a single rank.
* comments
* Container: Add `const_iterator`
* Fix MPI Comms
* Cleaning
* Remove PrintAll Leftover
* 1-D Output vector
* Reduced Diags: Support LoadBalance
* Cleaning of "Definitions ()"
* Updating inputs for testing Line
* data type specification
* IO
* Update inputs
* Update inputs
* error in header. Send to IO CPU
* Apply suggestions from code review
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* moved rank communication out of tile loop
* change m_data_vector. IO particle count
* Fixed input for rename. Gather particle number
* Apply suggestions from code review
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Gather
* Changed data output method to pushing values on vector
MPI Gather and Gatherv for data
Tell Evolve to run Load Balance
Tell InitData to run Load Balance
Define output method by printing valid particles
NOTE! Needs cleaning, commenting, removing some debugging tools
* Apply suggestions from code review
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Suggestions from review
* defensive programming on if / if else
added comments throughout
removed temporary debugging lines
* Documentation for line detector option
* Whitespaces
* MPI_Gather -> amrex::ParallelDescriptor::Gather
* ParallelDiscriptor, vectors at the end, no more 1990's malloc for capacity allocation
* whitespaces
* output optimized for CSV
* Python notebook for reference
* Input file for current test on CORI
* 2D plane functionality
* Regression test
* Delete DoubleSlit_2021_11_17.ipynb
* Whitespace fix
* pandas
* Error set to 2.5%, fixed source
* style
* zenodo orcid
* Apply suggestions from code review
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
* Review changes and swapped MPI direct call for Amrex::ParallelDescriptor
* [pre-commit.ci] auto fixes from pre-commit.com hooks
for more information, see https://pre-commit.ci
* Update WarpX-tests.ini
Fix Regression Test
* Update WarpX-tests.ini
Open PMD in cmakeSetupOpts
* Update dependencies.rst
* Update analysis_field_probe.py
* Update WarpX-tests.ini
* Analysis Script: Executable
```
chmod a+x scriptname.py
```
and use explicitly `python3`
* openPMD: optional for this test
* Inputs: add `geometry.dims = 2`
* Remove: diag1.write_species = 0
- segfaults for plotfiles (bug?)
- not needed, since we have no particles anyway
* Fix: typo in analysis
* test requirements: pandas
* Fix: Types
* as string: `<red_diag>.probe_geometry`
change this to a string, which is more user-friendly
* Python Script: Simplify + Style
* C++: Clean Up
* Azure: Run `apt update`
Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>
Co-authored-by: David Grote <dpgrote@lbl.gov>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Diffstat (limited to 'Source')
-rw-r--r-- | Source/Diagnostics/ReducedDiags/FieldProbe.H | 40 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/FieldProbe.cpp | 364 | ||||
-rw-r--r-- | Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp | 3 | ||||
-rw-r--r-- | Source/Evolve/WarpXEvolve.cpp | 1 |
4 files changed, 315 insertions, 93 deletions
diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.H b/Source/Diagnostics/ReducedDiags/FieldProbe.H index 8028d987d..15c5d18f2 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.H +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.H @@ -17,6 +17,17 @@ #include <unordered_map> #include <string> +#include <vector> + +/** + * This enumeration is used for assigning structural geometry levels (point vs line vs plane) + */ +enum struct DetectorGeometry +{ + Point = 0, + Line, + Plane +}; /** * This class mainly contains a function that computes the value of each component @@ -52,18 +63,41 @@ public: * Define constants used throughout FieldProbe */ - //! noutputs is 7 (Ex, Ey, Ez, Bx, By, Bz, S) - static constexpr int noutputs = FieldProbePIdx::nattribs; + //! noutputs is 10 (x, y, z, Ex, Ey, Ez, Bx, By, Bz, S) + static constexpr int noutputs = FieldProbePIdx::nattribs + 3; private: amrex::Real x_probe, y_probe, z_probe; + amrex::Real x1_probe, y1_probe, z1_probe; + amrex::Real target_normal_x, target_normal_y, target_normal_z; + amrex::Real target_up_x, target_up_y, target_up_z; + amrex::Real detector_radius; + + //! counts number of particles for all MPI ranks + long m_valid_particles {0}; + + //! determines geometry of detector point distribution + DetectorGeometry m_probe_geometry = DetectorGeometry::Point; + + //! determines number of particles places for non-point geometries + int m_resolution = 0; + + //! Empty vector for to which data is pushed + amrex::Vector<amrex::Real> m_data; + + //! Empty array to be used by IOProcessor node to store and output data + amrex::Vector<amrex::Real> m_data_out; //! this is the particle container in which probe particles are stored FieldProbeParticleContainer m_probe; - //! if true, integrate values over time instead of probing instantaneous values + //! if true, integrate values over time instead of probing instantaneous values bool m_field_probe_integrate = false; + + //! particle shape used for field gather int interp_order = 1; + + //! Judges whether to gather raw fields or interpolated data bool raw_fields = false; /** diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index d29795c35..3abf83347 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -31,9 +31,13 @@ #include <AMReX_Vector.H> #include <cmath> +#include <cstddef> +#include <cstdlib> +#include <iostream> #include <ostream> #include <string> #include <unordered_map> +#include <vector> using namespace amrex; @@ -72,18 +76,55 @@ FieldProbe::FieldProbe (std::string rd_name) * Define whether ot not to integrate fields */ amrex::ParmParse pp_rd_name(rd_name); - getWithParser(pp_rd_name, "x_probe", x_probe); - getWithParser(pp_rd_name, "y_probe", y_probe); - getWithParser(pp_rd_name, "z_probe", z_probe); + std::string m_probe_geometry_str = "Point"; + pp_rd_name.query("probe_geometry", m_probe_geometry_str); + + if (m_probe_geometry_str == "Point") + { + m_probe_geometry = DetectorGeometry::Point; + getWithParser(pp_rd_name, "x_probe", x_probe); + getWithParser(pp_rd_name, "y_probe", y_probe); + getWithParser(pp_rd_name, "z_probe", z_probe); + } + else if (m_probe_geometry_str == "Line") + { + m_probe_geometry = DetectorGeometry::Line; + getWithParser(pp_rd_name, "x_probe", x_probe); + getWithParser(pp_rd_name, "y_probe", y_probe); + getWithParser(pp_rd_name, "z_probe", z_probe); + getWithParser(pp_rd_name, "x1_probe", x1_probe); + getWithParser(pp_rd_name, "y1_probe", y1_probe); + getWithParser(pp_rd_name, "z1_probe", z1_probe); + getWithParser(pp_rd_name, "resolution", m_resolution); + } + else if (m_probe_geometry_str == "Plane") + { + m_probe_geometry = DetectorGeometry::Plane; + getWithParser(pp_rd_name, "x_probe", x_probe); + getWithParser(pp_rd_name, "y_probe", y_probe); + getWithParser(pp_rd_name, "z_probe", z_probe); + getWithParser(pp_rd_name, "target_normal_x", target_normal_x); + getWithParser(pp_rd_name, "target_normal_y", target_normal_y); + getWithParser(pp_rd_name, "target_normal_z", target_normal_z); + getWithParser(pp_rd_name, "target_up_x", target_up_x); + getWithParser(pp_rd_name, "target_up_y", target_up_y); + getWithParser(pp_rd_name, "target_up_z", target_up_z); + getWithParser(pp_rd_name, "detector_radius", detector_radius); + getWithParser(pp_rd_name, "resolution", m_resolution); + } + else + { + std::string err_str = "ERROR: Invalid probe geometry '"; + err_str.append(m_probe_geometry_str); + err_str.append("'. Valid geometries are Point, Line or Plane."); + amrex::Abort(err_str); + } pp_rd_name.query("integrate", m_field_probe_integrate); pp_rd_name.query("raw_fields", raw_fields); pp_rd_name.query("interp_order", interp_order); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(interp_order <= WarpX::nox , "Field probe interp_order should be less than or equal to algo.particle_shape"); - // resize data array - m_data.resize(noutputs, 0.0_rt); - if (ParallelDescriptor::IOProcessor()) { if ( m_IsNotRestart ) @@ -93,7 +134,6 @@ FieldProbe::FieldProbe (std::string rd_name) // write header row int c = 0; - ofs << "#"; ofs << "[" << c++ << "]step()"; ofs << m_sep; ofs << "[" << c++ << "]time(s)"; @@ -104,44 +144,50 @@ FieldProbe::FieldProbe (std::string rd_name) { u_map = { - {FieldProbePIdx::Ex , " (V*s/m) "}, - {FieldProbePIdx::Ey , " (V*s/m) "}, - {FieldProbePIdx::Ez , " (V*s/m) "}, - {FieldProbePIdx::Bx , " (T*s) "}, - {FieldProbePIdx::By , " (T*s) "}, - {FieldProbePIdx::Bz , " (T*s) "}, - {FieldProbePIdx::S , " (W*s/m^2) "} + {FieldProbePIdx::Ex , "-(V*s/m)"}, + {FieldProbePIdx::Ey , "-(V*s/m)"}, + {FieldProbePIdx::Ez , "-(V*s/m)"}, + {FieldProbePIdx::Bx , "-(T*s)"}, + {FieldProbePIdx::By , "-(T*s)"}, + {FieldProbePIdx::Bz , "-(T*s)"}, + {FieldProbePIdx::S , "-(W*s/m^2)"} }; } else { u_map = { - {FieldProbePIdx::Ex , " (V/m) "}, - {FieldProbePIdx::Ey , " (V/m) "}, - {FieldProbePIdx::Ez , " (V/m) "}, - {FieldProbePIdx::Bx , " (T) "}, - {FieldProbePIdx::By , " (T) "}, - {FieldProbePIdx::Bz , " (T) "}, - {FieldProbePIdx::S , " (W/m^2) "} + {FieldProbePIdx::Ex , "-(V/m)"}, + {FieldProbePIdx::Ey , "-(V/m)"}, + {FieldProbePIdx::Ez , "-(V/m)"}, + {FieldProbePIdx::Bx , "-(T)"}, + {FieldProbePIdx::By , "-(T)"}, + {FieldProbePIdx::Bz , "-(T)"}, + {FieldProbePIdx::S , "-(W/m^2)"} }; } for (int lev = 0; lev < nLevel; ++lev) { ofs << m_sep; - ofs << "[" << c++ << "]probe_Ex_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Ex]; + ofs << "[" << c++ << "]part_x_lev" + std::to_string(lev) + "-(m)"; + ofs << m_sep; + ofs << "[" << c++ << "]part_y_lev" + std::to_string(lev) + "-(m)"; + ofs << m_sep; + ofs << "[" << c++ << "]part_z_lev" + std::to_string(lev) + "-(m)"; + ofs << m_sep; + ofs << "[" << c++ << "]part_Ex_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Ex]; ofs << m_sep; - ofs << "[" << c++ << "]probe_Ey_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Ey]; + ofs << "[" << c++ << "]part_Ey_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Ey]; ofs << m_sep; - ofs << "[" << c++ << "]probe_Ez_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Ez]; + ofs << "[" << c++ << "]part_Ez_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Ez]; ofs << m_sep; - ofs << "[" << c++ << "]probe_Bx_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Bx]; + ofs << "[" << c++ << "]part_Bx_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Bx]; ofs << m_sep; - ofs << "[" << c++ << "]probe_By_lev" + std::to_string(lev) + u_map[FieldProbePIdx::By]; + ofs << "[" << c++ << "]part_By_lev" + std::to_string(lev) + u_map[FieldProbePIdx::By]; ofs << m_sep; - ofs << "[" << c++ << "]probe_Bz_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Bz]; + ofs << "[" << c++ << "]part_Bz_lev" + std::to_string(lev) + u_map[FieldProbePIdx::Bz]; ofs << m_sep; - ofs << "[" << c++ << "]probe_S_lev" + std::to_string(lev) + u_map[FieldProbePIdx::S]; //update all units if integrating (might be energy / m^2) + ofs << "[" << c++ << "]part_S_lev" + std::to_string(lev) + u_map[FieldProbePIdx::S]; } ofs << std::endl; @@ -153,24 +199,111 @@ FieldProbe::FieldProbe (std::string rd_name) void FieldProbe::InitData () { - //create 1D array for X, Y, and Z of particles - amrex::Vector<amrex::ParticleReal> xpos; - amrex::Vector<amrex::ParticleReal> ypos; - amrex::Vector<amrex::ParticleReal> zpos; + if (m_probe_geometry == DetectorGeometry::Point) + { - // for now, only one MPI rank adds a probe particle - if (ParallelDescriptor::IOProcessor()) + // create 1D vector for X, Y, and Z of particles + amrex::Vector<amrex::ParticleReal> xpos; + amrex::Vector<amrex::ParticleReal> ypos; + amrex::Vector<amrex::ParticleReal> zpos; + + // for now, only one MPI rank adds a probe particle + if (ParallelDescriptor::IOProcessor()) + { + xpos.push_back(x_probe); + ypos.push_back(y_probe); + zpos.push_back(z_probe); + } + + // add particles on lev 0 to m_probe + m_probe.AddNParticles(0, xpos, ypos, zpos); + } + else if (m_probe_geometry == DetectorGeometry::Line) { - xpos.push_back(x_probe); - ypos.push_back(y_probe); - zpos.push_back(z_probe); + amrex::Vector<amrex::ParticleReal> xpos; + amrex::Vector<amrex::ParticleReal> ypos; + amrex::Vector<amrex::ParticleReal> zpos; + if (ParallelDescriptor::IOProcessor()) + { + xpos.reserve(m_resolution); + ypos.reserve(m_resolution); + zpos.reserve(m_resolution); + + // Final - initial / steps. Array contains dx, dy, dz + amrex::Real DetLineStepSize[3]{ + (x1_probe - x_probe) / (m_resolution - 1), + (y1_probe - y_probe) / (m_resolution - 1), + (z1_probe - z_probe) / (m_resolution - 1)}; + for ( int step = 0; step < m_resolution; step++) + { + xpos.push_back(x_probe + (DetLineStepSize[0] * step)); + ypos.push_back(y_probe + (DetLineStepSize[1] * step)); + zpos.push_back(z_probe + (DetLineStepSize[2] * step)); + } + } + m_probe.AddNParticles(0, xpos, ypos, zpos); + } + else if (m_probe_geometry == DetectorGeometry::Plane) + { + amrex::Vector<amrex::ParticleReal> xpos; + amrex::Vector<amrex::ParticleReal> ypos; + amrex::Vector<amrex::ParticleReal> zpos; + if (ParallelDescriptor::IOProcessor()) + { + std::size_t const res2 = std::size_t(m_resolution) * std::size_t(m_resolution); + xpos.reserve(res2); + ypos.reserve(res2); + zpos.reserve(res2); + + // create vector orthonormal to input vectors + amrex::Real orthotarget[3]{ + target_normal_y * target_up_z - target_normal_z * target_up_y, + target_normal_z * target_up_x - target_normal_x * target_up_z, + target_normal_x * target_up_y - target_normal_y * target_up_x}; + // find upper left and lower right bounds of detector + amrex::Real direction[3]{ + orthotarget[0] - target_up_x, + orthotarget[1] - target_up_y, + orthotarget[2] - target_up_z}; + amrex::Real upperleft[3]{ + x_probe - (direction[0] * detector_radius), + y_probe - (direction[1] * detector_radius), + z_probe - (direction[2] * detector_radius)}; + amrex::Real lowerright[3]{ + x_probe + (direction[0] * detector_radius), + y_probe + (direction[1] * detector_radius), + z_probe + (direction[2] * detector_radius)}; + // create array containing point-to-point step size + amrex::Real DetPlaneStepSize[3]{ + (lowerright[0] - upperleft[0]) / (m_resolution - 1), + (lowerright[1] - upperleft[1]) / (m_resolution - 1), + (lowerright[2] - upperleft[2]) / (m_resolution - 1)}; + amrex::Real temp_pos[3]{}; + // Target point on top of plane (arbitrarily top of plane perpendicular to yz) + // For each point along top of plane, fill in YZ's beneath, then push back + for ( int step = 0; step < m_resolution; step++) + { + temp_pos[0] = upperleft[0] + (DetPlaneStepSize[0] * step); + for ( int yzstep = 0; yzstep < m_resolution; yzstep++) + { + temp_pos[1] = upperleft[1] + (DetPlaneStepSize[1] * yzstep); + temp_pos[2] = upperleft[2] + (DetPlaneStepSize[2] * yzstep); + xpos.push_back(temp_pos[0]); + ypos.push_back(temp_pos[1]); + zpos.push_back(temp_pos[2]); + } + } + } + m_probe.AddNParticles(0, xpos, ypos, zpos); + } + else + { + amrex::Abort("ERROR: Invalid probe geometry. Valid geometries are point (0) and line (1)."); } - - // add np particles on lev 0 to m_probe - m_probe.AddNParticles(0, xpos, ypos, zpos); } -void FieldProbe::LoadBalance() { +void FieldProbe::LoadBalance () +{ m_probe.Redistribute(); } @@ -239,15 +372,28 @@ void FieldProbe::ComputeDiags (int step) amrex::IndexType const Bytype = By.ixType(); amrex::IndexType const Bztype = Bz.ixType(); - //defined for use in determining which CPU contains the particles - int probe_proc = -1; - // loop over each particle // TODO: add OMP parallel as in PhysicalParticleContainer::Evolve + long numparticles = 0; // particles on this MPI rank using MyParIter = FieldProbeParticleContainer::iterator; for (MyParIter pti(m_probe, lev); pti.isValid(); ++pti) { + // count particle on MPI rank + numparticles += pti.numParticles(); + } + + if (m_intervals.contains(step+1)) + { + // reset m_data vector to clear pushed values. Reserves data + m_data.clear(); + m_data.shrink_to_fit(); + m_data.reserve(numparticles * noutputs); + } + + for (MyParIter pti(m_probe, lev); pti.isValid(); ++pti) + { const auto getPosition = GetParticlePosition(pti); + auto const np = pti.numParticles(); if( ProbeInDomain() ) { @@ -295,7 +441,6 @@ void FieldProbe::ComputeDiags (int step) const amrex::GpuArray<amrex::Real, 3> xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; const Dim3 lo = lbound(box); - // Interpolating to the probe positions for each particle // Temporarily defining modes and interp outside ParallelFor to avoid GPU compilation errors. const int temp_modes = WarpX::n_rz_azimuthal_modes; const int temp_interp_order = interp_order; @@ -303,7 +448,7 @@ void FieldProbe::ComputeDiags (int step) const bool temp_field_probe_integrate = m_field_probe_integrate; amrex::Real const dt = WarpX::GetInstance().getdt(lev); - long const np = pti.numParticles(); + // Interpolating to the probe positions for each particle amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long ip) { amrex::ParticleReal xp, yp, zp; @@ -354,7 +499,6 @@ void FieldProbe::ComputeDiags (int step) } else { - // Either save the interpolated fields or the raw fields depending on the raw_fields flag part_Ex[ip] = Exp; //remember to add lorentz transform part_Ey[ip] = Eyp; //remember to add lorentz transform part_Ez[ip] = Ezp; //remember to add lorentz transform @@ -364,63 +508,107 @@ void FieldProbe::ComputeDiags (int step) part_S[ip] = S; //remember to add lorentz transform } });// ParallelFor Close - // this check is here because for m_field_probe_integrate == True, we always compute // but we only write when we truly are in an output interval step - if (m_intervals.contains(step+1)) { - for (int ip = 0; ip < np; ip++) { - // Fill output array - m_data[ip * noutputs + FieldProbePIdx::Ex] = part_Ex[ip]; - m_data[ip * noutputs + FieldProbePIdx::Ey] = part_Ey[ip]; - m_data[ip * noutputs + FieldProbePIdx::Ez] = part_Ez[ip]; - m_data[ip * noutputs + FieldProbePIdx::Bx] = part_Bx[ip]; - m_data[ip * noutputs + FieldProbePIdx::By] = part_By[ip]; - m_data[ip * noutputs + FieldProbePIdx::Bz] = part_Bz[ip]; - m_data[ip * noutputs + FieldProbePIdx::S] = part_S[ip]; + if (m_intervals.contains(step+1)) + { + for (auto ip=0; ip < np; ip++) + { + amrex::ParticleReal xp, yp, zp; + getPosition(ip, xp, yp, zp); + + // push to output vector + m_data.push_back(xp); + m_data.push_back(yp); + m_data.push_back(zp); + m_data.push_back(part_Ex[ip]); + m_data.push_back(part_Ey[ip]); + m_data.push_back(part_Ez[ip]); + m_data.push_back(part_Bx[ip]); + m_data.push_back(part_By[ip]); + m_data.push_back(part_Bz[ip]); + m_data.push_back(part_S[ip]); } - /* m_data now contains up-to-date values for: - * [Ex, Ey, Ez, Bx, By, Bz, and S] */ + /* m_data now contains up-to-date values for: + * [x, y, z, Ex, Ey, Ez, Bx, By, Bz, and S] */ } - - // do we have the one and only probe particle on our MPI rank? - if (np > 0) - probe_proc = amrex::ParallelDescriptor::MyProc(); } - } // end particle iterator loop - - // make sure data is in m_data Gpu::synchronize(); + if (m_intervals.contains(step+1)) + { + // returns total number of mpi notes into mpisize + int mpisize = ParallelDescriptor::NProcs(); - // this check is here because for m_field_probe_integrate == True, we always compute - // but we only write when we truly are in an output interval step - if (m_intervals.contains(step+1)) { - /* - * All the processors have probe_proc = -1 except the one that contains the point, which - * has probe_proc equal to a number >=0. Therefore, ReduceIntMax communicates to all the - * processors the rank of the processor which contains the point - */ - amrex::ParallelDescriptor::ReduceIntMax(probe_proc); - - if (probe_proc != amrex::ParallelDescriptor::IOProcessorNumber() and probe_proc != -1) { - if (amrex::ParallelDescriptor::MyProc() == probe_proc) { - amrex::ParallelDescriptor::Send(m_data.data(), noutputs, - amrex::ParallelDescriptor::IOProcessorNumber(), - 0); - } - if (amrex::ParallelDescriptor::MyProc() - == amrex::ParallelDescriptor::IOProcessorNumber()) { - amrex::ParallelDescriptor::Recv(m_data.data(), noutputs, probe_proc, 0); + // allocates data space for length_array. Will contain size of m_data from each processor + amrex::Vector<int> length_vector; + amrex::Vector<int> localsize; + + if (amrex::ParallelDescriptor::IOProcessor()) { + length_vector.resize(mpisize, 0); + } + localsize.resize(1,0); + localsize[0] = m_data.size(); + + // gather size of m_data from each processor + amrex::ParallelDescriptor::Gather(localsize.data(), 1, + length_vector.data(), 1, + amrex::ParallelDescriptor::IOProcessorNumber()); + + // IO processor sums values from length_array to get size of total output array. + /* displs records the size of each m_data as well as previous displs. This array + * tells Gatherv where in the m_data_out array allocation to write incomming data. */ + long total_data_size = 0; + amrex::Vector<int> displs_vector; + if (amrex::ParallelDescriptor::IOProcessor()) { + displs_vector.resize(mpisize, 0); + displs_vector[0] = 0; + total_data_size += length_vector[0]; + for (int i=1; i<mpisize; i++) { + displs_vector[i] = (displs_vector[i-1] + length_vector[i-1]); + total_data_size += length_vector[i]; } + // valid particles are counted (for all MPI ranks) to inform output processes as to size of output + m_valid_particles = total_data_size / noutputs; + m_data_out.resize(total_data_size, 0); } - } // send to IO Processor + // resize receive buffer (resize, initialize 0) + // gather m_data of varied lengths from all processors. Prints to m_data_out + amrex::ParallelDescriptor::Gatherv(m_data.data(), localsize[0], + m_data_out.data(), length_vector, displs_vector, + amrex::ParallelDescriptor::IOProcessorNumber()); + } }// end loop over refinement levels + // make sure data is in m_data on the IOProcessor + // TODO: In the future, we want to use a parallel I/O method instead (plotfiles or openPMD) } // end void FieldProbe::ComputeDiags void FieldProbe::WriteToFile (int step) const { if (ProbeInDomain() && amrex::ParallelDescriptor::IOProcessor()) { - ReducedDiags::WriteToFile (step); + // open file + std::ofstream ofs{m_path + m_rd_name + "." + m_extension, + std::ofstream::out | std::ofstream::app}; + + // loop over num valid particles and write + for (int i = 0; i < m_valid_particles; i++) + { + ofs << std::fixed << std::defaultfloat; + ofs << step + 1; + ofs << m_sep; + ofs << std::fixed << std::setprecision(14) << std::scientific; + // write time + ofs << WarpX::GetInstance().gett_new(0); + + for (int k = 0; k < noutputs; k++) + { + ofs << m_sep; + ofs << m_data_out[i * noutputs + k]; + } + ofs << std::endl; + } // end loop over data size + // close file + ofs.close(); } } diff --git a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp index c72362b6b..111e97e96 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbeParticleContainer.cpp @@ -108,8 +108,7 @@ FieldProbeParticleContainer::AddNParticles (int lev, #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(y) ; p.pos(0) = x[i]; - p.pos(1) = 0; - p.pos(2) = z[i]; + p.pos(1) = z[i]; #endif // write position, cpu id, and particle id to particle pinned_tile.push_back(p); diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 547f632e9..c3c297b56 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -335,6 +335,7 @@ WarpX::Evolve (int numsteps) /// reduced diags if (reduced_diags->m_plot_rd != 0) { + reduced_diags->LoadBalance(); reduced_diags->ComputeDiags(step); reduced_diags->WriteToFile(step); } |