aboutsummaryrefslogtreecommitdiff
path: root/Source
diff options
context:
space:
mode:
Diffstat (limited to 'Source')
-rw-r--r--Source/Diagnostics/BoostedFrameDiagnostic.H282
-rw-r--r--Source/Diagnostics/BoostedFrameDiagnostic.cpp729
-rw-r--r--Source/Initialization/WarpXInitData.cpp11
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp47
-rw-r--r--Source/Particles/Pusher/Make.package1
-rw-r--r--Source/Particles/Pusher/UpdateMomentumHigueraCary.H59
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.cpp14
-rw-r--r--Source/Utils/WarpXAlgorithmSelection.H3
-rw-r--r--Source/Utils/WarpXAlgorithmSelection.cpp1
-rw-r--r--Source/WarpX.H4
-rw-r--r--Source/WarpX.cpp40
11 files changed, 896 insertions, 295 deletions
diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H
index 4263af8d0..5d95aaf7d 100644
--- a/Source/Diagnostics/BoostedFrameDiagnostic.H
+++ b/Source/Diagnostics/BoostedFrameDiagnostic.H
@@ -13,116 +13,231 @@
#include "MultiParticleContainer.H"
#include "WarpXConst.H"
-///
-/// BoostedFrameDiagnostic is for handling IO when running in a boosted
-/// frame of reference. Because of the relativity of simultaneity, events that
-/// are synchronized in the simulation frame are not synchronized in the
-/// lab frame. Thus, at a given t_boost, we must write slices of data to
-/// multiple output files, each one corresponding to a given time in the lab frame.
-///
-class BoostedFrameDiagnostic {
- ///
- /// LabSnapShot stores metadata corresponding to a single time
- /// snapshot in the lab frame. The snapshot is written to disk
- /// in the directory "file_name". zmin_lab, zmax_lab, and t_lab
- /// are all constant for a given snapshot. current_z_lab and
- /// current_z_boost for each snapshot are updated as the
- /// simulation time in the boosted frame advances.
- ///
- struct LabSnapShot {
-
- std::string file_name;
- amrex::Real t_lab;
- amrex::RealBox prob_domain_lab_;
- amrex::IntVect prob_ncells_lab_;
-
- amrex::Real current_z_lab;
- amrex::Real current_z_boost;
-
- int ncomp_to_dump_;
- std::vector<std::string> mesh_field_names_;
-
- int file_num;
- int initial_i;
- const BoostedFrameDiagnostic& my_bfd;
-
- LabSnapShot(amrex::Real t_lab_in, amrex::Real t_boost,
- const amrex::RealBox prob_domain_lab,
- const amrex::IntVect prob_ncells_lab,
- const int ncomp_to_dump,
- const std::vector<std::string> mesh_field_names,
- int file_num_in,
- const BoostedFrameDiagnostic& bfd);
-
- ///
- /// This snapshot is at time t_lab, and the simulation is at time t_boost.
- /// The Lorentz transformation picks out one slice corresponding to both
- /// of those times, at position current_z_boost and current_z_lab in the
- /// boosted and lab frames, respectively.
- ///
- void updateCurrentZPositions(amrex::Real t_boost, amrex::Real inv_gamma,
- amrex::Real inv_beta);
-
- ///
- /// Write some useful metadata about this snapshot.
- ///
- void writeSnapShotHeader();
- };
+/** \brief
+ * The capability for back-transformed lab-frame data is implemented to generate
+ * the full diagnostic snapshot for the entire domain and reduced diagnostic
+ * (1D, 2D or 3D 'slices') for a sub-domain.
+ * LabFrameDiag class defines the parameters required to backtrasform data from
+ * boosted frame at (z_boost,t_boost) to lab-frame at (z_lab, t_lab) using Lorentz
+ * transformation. This Lorentz transformation picks out one slice corresponding
+ * to both of those times, at position current_z_boost and current_z_lab in the
+ * boosted and lab frames, respectively.
+ * Two derived classes, namely, LabFrameSnapShot and LabFrameSlice are defined to
+ * store the full back-transformed diagnostic snapshot of the entire domain and
+ * reduced back-transformed diagnostic for a sub-domain, respectively.
+ * The approach here is to define an array of LabFrameDiag which would include
+ * both, full domain snapshots and reduced domain 'slices', sorted based on their
+ * respective t_lab. This is done to re-use the backtransformed data stored in
+ * the slice multifab at (z_lab,t_lab)
+ * for the full domain snapshot and sub-domain slices that have the same t_lab,
+ * instead of re-generating the backtransformed slice data at z_lab for a given
+ * t_lab for each diagnostic.
+ */
- amrex::Real gamma_boost_;
+class LabFrameDiag {
+ public:
+ std::string file_name;
+ amrex::Real t_lab;
+ amrex::RealBox prob_domain_lab_;
+ amrex::IntVect prob_ncells_lab_;
+ amrex::RealBox diag_domain_lab_;
+ amrex::Box buff_box_;
+
+ amrex::Real current_z_lab;
+ amrex::Real current_z_boost;
amrex::Real inv_gamma_boost_;
- amrex::Real beta_boost_;
amrex::Real inv_beta_boost_;
amrex::Real dz_lab_;
- amrex::Real inv_dz_lab_;
- amrex::Real dt_snapshots_lab_;
- amrex::Real dt_boost_;
- int N_snapshots_;
- int boost_direction_;
+ amrex::Real dx_;
+ amrex::Real dy_;
- // For back-transformed diagnostics of grid fields, data_buffer_[i]
+ int ncomp_to_dump_;
+ std::vector<std::string> mesh_field_names_;
+
+ int file_num;
+ int initial_i;
+
+ // For back-transformed diagnostics of grid fields, data_buffer_
// stores a buffer of the fields in the lab frame (in a MultiFab, i.e.
// with all box data etc.). When the buffer if full, dump to file.
- amrex::Vector<std::unique_ptr<amrex::MultiFab> > data_buffer_;
+ std::unique_ptr<amrex::MultiFab> data_buffer_;
// particles_buffer_ is currently blind to refinement level.
- // particles_buffer_[i][j] is a WarpXParticleContainer::DiagnosticParticleData where
- // - i is the back-transformed snapshot number
- // - j is the species number
- amrex::Vector<amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> > particles_buffer_;
+ // particles_buffer_[j] is a WarpXParticleContainer::DiagnosticParticleData
+ // where - j is the species number for the current diag
+ amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> particles_buffer_;
+ // buff_counter_ is the number of z slices in data_buffer_
+ int buff_counter_;
int num_buffer_ = 256;
- int max_box_size_ = 256;
- // buff_counter_[i] is the number of z slices in data_buffer_[i]
- // for snapshot number i.
- amrex::Vector<int> buff_counter_;
+ int max_box_size = 256;
+ void updateCurrentZPositions(amrex::Real t_boost, amrex::Real inv_gamma,
+ amrex::Real inv_beta);
- amrex::Vector<LabSnapShot> snapshots_;
+ void createLabFrameDirectories();
- void writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata,
- const std::string& name, const int i_lab);
+ void writeLabFrameHeader();
+
+ /// Back-transformed lab-frame field data is copied from
+ /// tmp_slice to data_buffer where it is stored.
+ /// For the full diagnostic, all the data in the MultiFab is copied.
+ /// For the reduced diagnostic, the data is selectively copied if the
+ /// extent of the z_lab multifab intersects with the user-defined sub-domain
+ /// for the reduced diagnostic (i.e., a 1D, 2D, or 3D region of the domain).
+ virtual void AddDataToBuffer(amrex::MultiFab& tmp_slice_ptr, int i_lab,
+ amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump){};
+
+ /// Back-transformed lab-frame particles is copied from
+ /// tmp_particle_buffer to particles_buffer.
+ /// For the full diagnostic, all the particles are copied,
+ /// while for the reduced diagnostic, particles are selectively
+ /// copied if their position in within the user-defined
+ /// sub-domain +/- 1 cell size width for the reduced slice diagnostic.
+ virtual void AddPartDataToParticleBuffer(
+ amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer,
+ int nSpeciesBoostedFrame) {};
+};
+
+/** \brief
+ * LabFrameSnapShot stores the back-transformed lab-frame metadata
+ * corresponding to a single time snapshot of the full domain.
+ * The snapshot data is written to disk in the directory lab_frame_data/snapshots/.
+ * zmin_lab, zmax_lab, and t_lab are all constant for a given snapshot.
+ * current_z_lab and current_z_boost for each snapshot are updated as the
+ * simulation time in the boosted frame advances.
+ */
+
+class LabFrameSnapShot : public LabFrameDiag {
+ public:
+ LabFrameSnapShot(amrex::Real t_lab_in, amrex::Real t_boost,
+ amrex::Real inv_gamma_boost_in, amrex::Real inv_beta_boost_in,
+ amrex::Real dz_lab_in, amrex::RealBox prob_domain_lab,
+ amrex::IntVect prob_ncells_lab, int ncomp_to_dump,
+ std::vector<std::string> mesh_field_names,
+ amrex::RealBox diag_domain_lab,
+ amrex::Box diag_box, int file_num_in);
+ void AddDataToBuffer( amrex::MultiFab& tmp_slice, int k_lab,
+ amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump) override;
+ void AddPartDataToParticleBuffer(
+ amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer,
+ int nSpeciesBoostedFrame) override;
+};
+
+
+/** \brief
+ * LabFrameSlice stores the back-transformed metadata corresponding to a single time at the
+ * user-defined slice location. This could be a 1D line, 2D slice, or 3D box
+ * (a reduced back-transformed diagnostic) within the computational
+ * domain, as defined in the input file by the user. The slice is written to disk in the
+ * lab_frame_data/slices. Similar to snapshots, zmin_lab, zmax_lab, and
+ * t_lab are constant for a given slice. current_z_lab and current_z_boost
+ * for each snapshot are updated as the sim time in boosted frame advances.
+ */
+class LabFrameSlice : public LabFrameDiag {
+ public:
+ LabFrameSlice(amrex::Real t_lab_in, amrex::Real t_boost,
+ amrex::Real inv_gamma_boost_in, amrex::Real inv_beta_boost_in,
+ amrex::Real dz_lab_in, amrex::RealBox prob_domain_lab,
+ amrex::IntVect prob_ncells_lab, int ncomp_to_dump,
+ std::vector<std::string> mesh_field_names,
+ amrex::RealBox diag_domain_lab,
+ amrex::Box diag_box, int file_num_in,
+ amrex::Real cell_dx, amrex::Real cell_dy);
+ void AddDataToBuffer( amrex::MultiFab& tmp_slice_ptr, int i_lab,
+ amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump) override;
+ void AddPartDataToParticleBuffer(
+ amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer,
+ int nSpeciesBoostedFrame) override;
+};
+
+/** \brief
+ * BoostedFrameDiagnostic class handles the back-transformation of data when
+ * running simulations in a boosted frame of reference to the lab-frame.
+ * Because of the relativity of simultaneity, events that are synchronized
+ * in the simulation boosted frame are not
+ * synchronized in the lab frame. Thus, at a given t_boost, we must write
+ * slices of back-transformed data to multiple output files, each one
+ * corresponding to a given time in the lab frame. The member function
+ * writeLabFrameData() orchestrates the operations required to
+ * Lorentz-transform data from boosted-frame to lab-frame and store them
+ * in the LabFrameDiag class, which writes out the field and particle data
+ * to the output directory. The functions Flush() and writeLabFrameData()
+ * are called at the end of the simulation and when the
+ * the buffer for data storage is full, respectively. The particle data
+ * is collected and written only if particle.do_boosted_frame_diagnostic = 1.
+ */
+class BoostedFrameDiagnostic {
-#ifdef WARPX_USE_HDF5
- void writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdata,
- const std::string& name, const std::string& species_name);
-#endif
public:
BoostedFrameDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab,
amrex::Real v_window_lab, amrex::Real dt_snapshots_lab,
- int N_snapshots, amrex::Real gamma_boost,
- amrex::Real t_boost, amrex::Real dt_boost, int boost_direction,
- const amrex::Geometry& geom);
+ int N_snapshots, amrex::Real dt_slice_snapshots_lab,
+ int N_slice_snapshots, amrex::Real gamma_boost,
+ amrex::Real t_boost, amrex::Real dt_boost,
+ int boost_direction, const amrex::Geometry& geom,
+ amrex::RealBox& slice_realbox);
+ /// Flush() is called at the end of the simulation when the buffers that contain
+ /// back-transformed lab-frame data even if they are not full.
void Flush(const amrex::Geometry& geom);
+ /** \brief
+ * The order of operations performed in writeLabFrameData is as follows :
+ * 1. Loops over the sorted back-transformed diags and for each diag
+ * steps 2-7 are performed
+ * 2. Based on t_lab and t_boost, obtain z_lab and z_boost.
+ * 3. Define data_buffer multifab that will store the data in the BT diag.
+ * 4. Define slice multifab at z_index that corresponds to z_boost and
+ * getslicedata using cell-centered data at z_index and its distribution map.
+ * 5. Lorentz transform data stored in slice from z_boost,t_Boost to z_lab,t_lab
+ * and store in slice multifab.
+ * 6. Generate a temporary slice multifab with distribution map of lab-frame
+ * data but at z_boost and
+ * ParallelCopy data from the slice multifab to the temporary slice.
+ * 7. Finally, AddDataToBuffer is called where the data from temporary slice
+ * is simply copied from tmp_slice(i,j,k_boost) to
+ * LabFrameDiagSnapshot(i,j,k_lab) for full BT lab-frame diagnostic
+ * OR from tmp_slice(i,j,k_boost) to
+ * LabFrameDiagSlice(i,j,k_lab) for the reduced slice diagnostic
+ * 8. Similarly, particles that crossed the z_boost plane are selected
+ * and lorentz-transformed to the lab-frame and copied to the full
+ * and reduce diagnostic and stored in particle_buffer.
+ */
void writeLabFrameData(const amrex::MultiFab* cell_centered_data,
const MultiParticleContainer& mypc,
const amrex::Geometry& geom,
const amrex::Real t_boost, const amrex::Real dt);
-
+ /// The metadata containg information on t_boost, num_snapshots, and Lorentz parameters.
void writeMetaData();
private:
+ amrex::Real gamma_boost_;
+ amrex::Real inv_gamma_boost_;
+ amrex::Real beta_boost_;
+ amrex::Real inv_beta_boost_;
+ amrex::Real dz_lab_;
+ amrex::Real inv_dz_lab_;
+ amrex::Real dt_snapshots_lab_;
+ amrex::Real dt_boost_;
+ int N_snapshots_;
+ int boost_direction_;
+ int N_slice_snapshots_;
+ amrex::Real dt_slice_snapshots_lab_;
+
+ int num_buffer_ = 256;
+ int max_box_size_ = 256;
+
+ std::vector<std::unique_ptr<LabFrameDiag> > LabFrameDiags_;
+
+ void writeParticleData(
+ const WarpXParticleContainer::DiagnosticParticleData& pdata,
+ const std::string& name, const int i_lab);
+
+#ifdef WARPX_USE_HDF5
+ void writeParticleDataHDF5(
+ const WarpXParticleContainer::DiagnosticParticleData& pdata,
+ const std::string& name, const std::string& species_name);
+#endif
// Map field names and component number in cell_centered_data
std::map<std::string, int> possible_fields_to_dump = {
{"Ex" , 0},
@@ -147,6 +262,7 @@ private:
"jx", "jy", "jz", "rho"};
int ncomp_to_dump = 10;
+
};
#endif
diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp
index 45343a0cb..297b4f5be 100644
--- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp
+++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp
@@ -1,6 +1,8 @@
#include <AMReX_MultiFabUtil.H>
+#include <AMReX_MultiFabUtil_C.H>
#include "BoostedFrameDiagnostic.H"
+#include "SliceDiagnostic.H"
#include "WarpX_f.H"
#include "WarpX.H"
@@ -16,7 +18,9 @@ using namespace amrex;
*/
namespace
{
- const std::vector<std::string> particle_field_names = {"w", "x", "y", "z", "ux", "uy", "uz"};
+
+ const std::vector<std::string> particle_field_names = {"w", "x", "y",
+ "z", "ux", "uy", "uz"};
/*
Creates the HDF5 file in truncate mode and closes it.
@@ -336,7 +340,8 @@ namespace
*/
void output_write_field(const std::string& file_path,
const std::string& field_path,
- const MultiFab& mf, const int comp)
+ const MultiFab& mf, const int comp,
+ const int lo_x, const int lo_y, const int lo_z)
{
BL_PROFILE("output_write_field");
@@ -374,8 +379,15 @@ namespace
// slab lo index and shape.
#if (AMREX_SPACEDIM == 3)
hsize_t slab_offsets[3], slab_dims[3];
+ int shift[3];
+ shift[0] = lo_x;
+ shift[1] = lo_y;
+ shift[2] = lo_z;
#else
hsize_t slab_offsets[2], slab_dims[2];
+ int shift[2];
+ shift[0] = lo_x;
+ shift[1] = lo_z;
#endif
hid_t slab_dataspace;
@@ -396,7 +408,7 @@ namespace
{
AMREX_ASSERT(lo_vec[idim] >= 0);
AMREX_ASSERT(hi_vec[idim] > lo_vec[idim]);
- slab_offsets[idim] = lo_vec[idim];
+ slab_offsets[idim] = lo_vec[idim] - shift[idim];
slab_dims[idim] = hi_vec[idim] - lo_vec[idim] + 1;
}
@@ -444,39 +456,14 @@ namespace
}
#endif
-namespace
+bool compare_tlab_uptr(const std::unique_ptr<LabFrameDiag>&a,
+ const std::unique_ptr<LabFrameDiag>&b)
{
- void
- CopySlice(MultiFab& tmp, MultiFab& buf, int k_lab,
- const Gpu::ManagedDeviceVector<int>& map_actual_fields_to_dump)
- {
- const int ncomp_to_dump = map_actual_fields_to_dump.size();
- // Copy data from MultiFab tmp to MultiFab data_buffer[i].
-#ifdef _OPENMP
-#pragma omp parallel if (Gpu::notInLaunchRegion())
-#endif
- for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
- Array4< Real> tmp_arr = tmp[mfi].array();
- Array4< Real> buf_arr = buf[mfi].array();
- // For 3D runs, tmp is a 2D (x,y) multifab, that contains only
- // slice to write to file.
- const Box& bx = mfi.tilebox();
-
- const auto field_map_ptr = map_actual_fields_to_dump.dataPtr();
- ParallelFor(bx, ncomp_to_dump,
- [=] AMREX_GPU_DEVICE (int i, int j, int k, int n)
- {
- const int icomp = field_map_ptr[n];
-#if (AMREX_SPACEDIM == 3)
- buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp);
-#else
- buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp);
-#endif
- }
- );
- }
- }
+ return a->t_lab < b->t_lab;
+}
+namespace
+{
void
LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp)
{
@@ -497,21 +484,27 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp)
// Transform the transverse E and B fields. Note that ez and bz are not
// changed by the tranform.
Real e_lab, b_lab, j_lab, r_lab;
- e_lab = gamma_boost * (arr(i, j, k, 0) + beta_boost*clight*arr(i, j, k, 4));
- b_lab = gamma_boost * (arr(i, j, k, 4) + beta_boost*arr(i, j, k, 0)/clight);
+ e_lab = gamma_boost * (arr(i, j, k, 0) +
+ beta_boost*clight*arr(i, j, k, 4));
+ b_lab = gamma_boost * (arr(i, j, k, 4) +
+ beta_boost*arr(i, j, k, 0)/clight);
arr(i, j, k, 0) = e_lab;
arr(i, j, k, 4) = b_lab;
- e_lab = gamma_boost * (arr(i, j, k, 1) - beta_boost*clight*arr(i, j, k, 3));
- b_lab = gamma_boost * (arr(i, j, k, 3) - beta_boost*arr(i, j, k, 1)/clight);
+ e_lab = gamma_boost * (arr(i, j, k, 1) -
+ beta_boost*clight*arr(i, j, k, 3));
+ b_lab = gamma_boost * (arr(i, j, k, 3) -
+ beta_boost*arr(i, j, k, 1)/clight);
arr(i, j, k, 1) = e_lab;
arr(i, j, k, 3) = b_lab;
// Transform the charge and current density. Only the z component of j is affected.
- j_lab = gamma_boost*(arr(i, j, k, 8) + beta_boost*clight*arr(i, j, k, 9));
- r_lab = gamma_boost*(arr(i, j, k, 9) + beta_boost*arr(i, j, k, 8)/clight);
+ j_lab = gamma_boost*(arr(i, j, k, 8) +
+ beta_boost*clight*arr(i, j, k, 9));
+ r_lab = gamma_boost*(arr(i, j, k, 9) +
+ beta_boost*arr(i, j, k, 8)/clight);
arr(i, j, k, 8) = j_lab;
arr(i, j, k, 9) = r_lab;
@@ -524,15 +517,20 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp)
BoostedFrameDiagnostic::
BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab,
Real dt_snapshots_lab, int N_snapshots,
+ Real dt_slice_snapshots_lab, int N_slice_snapshots,
Real gamma_boost, Real t_boost, Real dt_boost,
- int boost_direction, const Geometry& geom)
+ int boost_direction, const Geometry& geom,
+ amrex::RealBox& slice_realbox)
: gamma_boost_(gamma_boost),
dt_snapshots_lab_(dt_snapshots_lab),
dt_boost_(dt_boost),
N_snapshots_(N_snapshots),
- boost_direction_(boost_direction)
+ boost_direction_(boost_direction),
+ N_slice_snapshots_(N_slice_snapshots),
+ dt_slice_snapshots_lab_(dt_slice_snapshots_lab)
{
+
BL_PROFILE("BoostedFrameDiagnostic::BoostedFrameDiagnostic");
AMREX_ALWAYS_ASSERT(WarpX::do_boosted_frame_fields or
@@ -553,12 +551,8 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab,
// Ny_lab = 1;
IntVect prob_ncells_lab = {Nx_lab, Nz_lab};
#endif
-
writeMetaData();
- if (WarpX::do_boosted_frame_fields) data_buffer_.resize(N_snapshots);
- if (WarpX::do_boosted_frame_particles) particles_buffer_.resize(N_snapshots);
-
// Query fields to dump
std::vector<std::string> user_fields_to_dump;
ParmParse pp("warpx");
@@ -581,6 +575,9 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab,
}
}
+ // allocating array with total number of lab frame diags (snapshots+slices)
+ LabFrameDiags_.resize(N_snapshots+N_slice_snapshots);
+
for (int i = 0; i < N_snapshots; ++i) {
Real t_lab = i * dt_snapshots_lab_;
// Get simulation domain physical coordinates (in boosted frame).
@@ -589,14 +586,95 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab,
// x and y bounds are the same for lab frame and boosted frame
prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_lab);
prob_domain_lab.setHi(AMREX_SPACEDIM-1, zmax_lab + v_window_lab * t_lab);
- // Construct LabSnapShot
- LabSnapShot snapshot(t_lab, t_boost, prob_domain_lab,
- prob_ncells_lab, ncomp_to_dump,
- mesh_field_names, i, *this);
- snapshots_.push_back(snapshot);
- buff_counter_.push_back(0);
- if (WarpX::do_boosted_frame_fields) data_buffer_[i].reset( nullptr );
+ Box diag_box = geom.Domain();
+ LabFrameDiags_[i].reset(new LabFrameSnapShot(t_lab, t_boost,
+ inv_gamma_boost_, inv_beta_boost_, dz_lab_,
+ prob_domain_lab, prob_ncells_lab,
+ ncomp_to_dump, mesh_field_names, prob_domain_lab,
+ diag_box, i));
+ }
+
+
+ for (int i = 0; i < N_slice_snapshots; ++i) {
+
+ amrex::Real cell_dx = 0;
+ amrex::Real cell_dy = 0;
+ IntVect slice_ncells_lab ;
+
+ // To construct LabFrameSlice(), the location of lo() and hi() of the
+ // reduced diag is computed using the user-defined values of the
+ // reduced diag (1D, 2D, or 3D). For visualization of the diagnostics,
+ // the number of cells in each dimension is required and
+ // is computed below for the reduced back-transformed lab-frame diag,
+ // similar to the full-diag.
+ const amrex::Real* current_slice_lo = slice_realbox.lo();
+ const amrex::Real* current_slice_hi = slice_realbox.hi();
+
+ const amrex::Real zmin_slice_lab = current_slice_lo[AMREX_SPACEDIM-1] /
+ ( (1.+beta_boost_)*gamma_boost_);
+ const amrex::Real zmax_slice_lab = current_slice_hi[AMREX_SPACEDIM-1] /
+ ( (1.+beta_boost_)*gamma_boost_);
+ int Nz_slice_lab = static_cast<unsigned>((
+ zmax_slice_lab - zmin_slice_lab) * inv_dz_lab_);
+ int Nx_slice_lab = ( current_slice_hi[0] - current_slice_lo[0] ) /
+ geom.CellSize(0);
+ if (Nx_slice_lab == 0 ) Nx_slice_lab = 1;
+ // if the x-dimension is reduced, increase total_cells by 1
+ // to be consistent with the number of cells created for the output.
+ if (Nx_lab != Nx_slice_lab) Nx_slice_lab++;
+ cell_dx = geom.CellSize(0);
+#if (AMREX_SPACEDIM == 3)
+ int Ny_slice_lab = ( current_slice_hi[1] - current_slice_lo[1]) /
+ geom.CellSize(1);
+ if (Ny_slice_lab == 0 ) Ny_slice_lab = 1;
+ // if the y-dimension is reduced, increase total_cells by 1
+ // to be consistent with the number of cells created for the output.
+ if (Ny_lab != Ny_slice_lab) Ny_slice_lab++;
+ slice_ncells_lab = {Nx_slice_lab, Ny_slice_lab, Nz_slice_lab};
+ cell_dy = geom.CellSize(1);
+#else
+ slice_ncells_lab = {Nx_slice_lab, Nz_slice_lab};
+#endif
+
+ IntVect slice_lo(AMREX_D_DECL(0,0,0));
+ IntVect slice_hi(AMREX_D_DECL(1,1,1));
+
+ for ( int i_dim=0; i_dim<AMREX_SPACEDIM; ++i_dim)
+ {
+ slice_lo[i_dim] = (slice_realbox.lo(i_dim) - geom.ProbLo(i_dim) -
+ 0.5*geom.CellSize(i_dim))/geom.CellSize(i_dim);
+ slice_hi[i_dim] = (slice_realbox.hi(i_dim) - geom.ProbLo(i_dim) -
+ 0.5*geom.CellSize(i_dim))/geom.CellSize(i_dim);
+ if (slice_lo[i_dim] == slice_hi[i_dim])
+ {
+ slice_hi[i_dim] = slice_lo[i_dim] + 1;
+ }
+ }
+ Box stmp(slice_lo,slice_hi);
+ Box slicediag_box = stmp;
+
+ Real t_slice_lab = i * dt_slice_snapshots_lab_ ;
+ RealBox prob_domain_lab = geom.ProbDomain();
+ // replace z bounds by lab-frame coordinates
+ prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_slice_lab);
+ prob_domain_lab.setHi(AMREX_SPACEDIM-1, zmax_lab + v_window_lab * t_slice_lab);
+ RealBox slice_dom_lab = slice_realbox;
+ // replace z bounds of slice in lab-frame coordinates
+ // note : x and y bounds are the same for lab and boosted frames
+ // initial lab slice extent //
+ slice_dom_lab.setLo(AMREX_SPACEDIM-1, zmin_slice_lab + v_window_lab * t_slice_lab );
+ slice_dom_lab.setHi(AMREX_SPACEDIM-1, zmax_slice_lab +
+ v_window_lab * t_slice_lab );
+
+ // construct labframeslice
+ LabFrameDiags_[i+N_snapshots].reset(new LabFrameSlice(t_slice_lab, t_boost,
+ inv_gamma_boost_, inv_beta_boost_, dz_lab_,
+ prob_domain_lab, slice_ncells_lab,
+ ncomp_to_dump, mesh_field_names, slice_dom_lab,
+ slicediag_box, i, cell_dx, cell_dy));
}
+ // sort diags based on their respective t_lab
+ std::stable_sort(LabFrameDiags_.begin(), LabFrameDiags_.end(), compare_tlab_uptr);
AMREX_ALWAYS_ASSERT(max_box_size_ >= num_buffer_);
}
@@ -612,18 +690,19 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom)
const std::vector<std::string> species_names = mypc.GetSpeciesNames();
// Loop over BFD snapshots
- for (int i = 0; i < N_snapshots_; ++i) {
+ for (int i = 0; i < LabFrameDiags_.size(); ++i) {
- Real zmin_lab = snapshots_[i].prob_domain_lab_.lo(AMREX_SPACEDIM-1);
- int i_lab = (snapshots_[i].current_z_lab - zmin_lab) / dz_lab_;
+ Real zmin_lab = LabFrameDiags_[i]->prob_domain_lab_.lo(AMREX_SPACEDIM-1);
+ int i_lab = (LabFrameDiags_[i]->current_z_lab - zmin_lab) / dz_lab_;
- if (buff_counter_[i] != 0) {
+ if (LabFrameDiags_[i]->buff_counter_ != 0) {
if (WarpX::do_boosted_frame_fields) {
- const BoxArray& ba = data_buffer_[i]->boxArray();
+ const BoxArray& ba = LabFrameDiags_[i]->data_buffer_->boxArray();
const int hi = ba[0].bigEnd(boost_direction_);
- const int lo = hi - buff_counter_[i] + 1;
+ const int lo = hi - LabFrameDiags_[i]->buff_counter_ + 1;
- Box buff_box = geom.Domain();
+ //Box buff_box = geom.Domain();
+ Box buff_box = LabFrameDiags_[i]->buff_box_;
buff_box.setSmall(boost_direction_, lo);
buff_box.setBig(boost_direction_, hi);
@@ -631,18 +710,23 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom)
buff_ba.maxSize(max_box_size_);
DistributionMapping buff_dm(buff_ba);
- const int ncomp = data_buffer_[i]->nComp();
+ const int ncomp = LabFrameDiags_[i]->data_buffer_->nComp();
MultiFab tmp(buff_ba, buff_dm, ncomp, 0);
- tmp.copy(*data_buffer_[i], 0, 0, ncomp);
+ tmp.copy(*LabFrameDiags_[i]->data_buffer_, 0, 0, ncomp);
#ifdef WARPX_USE_HDF5
- for (int comp = 0; comp < ncomp; ++comp)
- output_write_field(snapshots_[i].file_name, mesh_field_names[comp], tmp, comp);
+ for (int comp = 0; comp < ncomp; ++comp) {
+ output_write_field(LabFrameDiags_[i]->file_name,
+ mesh_field_names[comp], tmp, comp,
+ lbound(buff_box).x, lbound(buff_box).y,
+ lbound(buff_box).z);
+ }
#else
std::stringstream ss;
- ss << snapshots_[i].file_name << "/Level_0/" << Concatenate("buffer", i_lab, 5);
+ ss << LabFrameDiags_[i]->file_name << "/Level_0/"
+ << Concatenate("buffer", i_lab, 5);
VisMF::Write(tmp, ss.str());
#endif
}
@@ -655,25 +739,31 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom)
species_names[mypc.mapSpeciesBoostedFrameDiags(j)];
#ifdef WARPX_USE_HDF5
// Dump species data
- writeParticleDataHDF5(particles_buffer_[i][j],
- snapshots_[i].file_name,
+ writeParticleDataHDF5(LabFrameDiags_[i]->particles_buffer_[j],
+ LabFrameDiags_[i]->file_name,
species_name);
#else
std::stringstream part_ss;
- part_ss << snapshots_[i].file_name + "/" + species_name + "/";
+ part_ss << LabFrameDiags_[i]->file_name + "/" +
+ species_name + "/";
// Dump species data
- writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab);
+ writeParticleData(LabFrameDiags_[i]->particles_buffer_[j],
+ part_ss.str(), i_lab);
#endif
}
- particles_buffer_[i].clear();
+ LabFrameDiags_[i]->particles_buffer_.clear();
}
- buff_counter_[i] = 0;
+ LabFrameDiags_[i]->buff_counter_ = 0;
}
}
VisMF::SetHeaderVersion(current_version);
}
+
+
+
+
void
BoostedFrameDiagnostic::
writeLabFrameData(const MultiFab* cell_centered_data,
@@ -681,7 +771,6 @@ writeLabFrameData(const MultiFab* cell_centered_data,
const Geometry& geom, const Real t_boost, const Real dt) {
BL_PROFILE("BoostedFrameDiagnostic::writeLabFrameData");
-
VisMF::Header::Version current_version = VisMF::GetHeaderVersion();
VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1);
@@ -690,97 +779,140 @@ writeLabFrameData(const MultiFab* cell_centered_data,
const Real zhi_boost = domain_z_boost.hi(boost_direction_);
const std::vector<std::string> species_names = mypc.GetSpeciesNames();
+ Real prev_t_lab = -dt;
+ std::unique_ptr<amrex::MultiFab> tmp_slice_ptr;
+ std::unique_ptr<amrex::MultiFab> slice;
+ amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer;
// Loop over snapshots
- for (int i = 0; i < N_snapshots_; ++i) {
-
+ for (int i = 0; i < LabFrameDiags_.size(); ++i) {
// Get updated z position of snapshot
- const Real old_z_boost = snapshots_[i].current_z_boost;
- snapshots_[i].updateCurrentZPositions(t_boost,
+ const Real old_z_boost = LabFrameDiags_[i]->current_z_boost;
+ LabFrameDiags_[i]->updateCurrentZPositions(t_boost,
inv_gamma_boost_,
inv_beta_boost_);
- Real zmin_lab = snapshots_[i].prob_domain_lab_.lo(AMREX_SPACEDIM-1);
- Real zmax_lab = snapshots_[i].prob_domain_lab_.hi(AMREX_SPACEDIM-1);
+ Real diag_zmin_lab = LabFrameDiags_[i]->diag_domain_lab_.lo(AMREX_SPACEDIM-1);
+ Real diag_zmax_lab = LabFrameDiags_[i]->diag_domain_lab_.hi(AMREX_SPACEDIM-1);
- // If snapshot out of the domain, nothing to do
- if ( (snapshots_[i].current_z_boost < zlo_boost) or
- (snapshots_[i].current_z_boost > zhi_boost) or
- (snapshots_[i].current_z_lab < zmin_lab) or
- (snapshots_[i].current_z_lab > zmax_lab) ) continue;
+ if ( ( LabFrameDiags_[i]->current_z_boost < zlo_boost) or
+ ( LabFrameDiags_[i]->current_z_boost > zhi_boost) or
+ ( LabFrameDiags_[i]->current_z_lab < diag_zmin_lab) or
+ ( LabFrameDiags_[i]->current_z_lab > diag_zmax_lab) ) continue;
// Get z index of data_buffer_ (i.e. in the lab frame) where
// simulation domain (t', [zmin',zmax']), back-transformed to lab
// frame, intersects with snapshot.
- int i_lab = (snapshots_[i].current_z_lab - zmin_lab) / dz_lab_;
-
+ Real dom_zmin_lab = LabFrameDiags_[i]->prob_domain_lab_.lo(AMREX_SPACEDIM-1);
+ int i_lab = ( LabFrameDiags_[i]->current_z_lab - dom_zmin_lab) / dz_lab_;
// If buffer of snapshot i is empty...
- if (buff_counter_[i] == 0) {
- // ... reset fields buffer data_buffer_[i]
+ if ( LabFrameDiags_[i]->buff_counter_ == 0) {
+ // ... reset fields buffer data_buffer_
if (WarpX::do_boosted_frame_fields) {
- Box buff_box = geom.Domain();
- buff_box.setSmall(boost_direction_, i_lab - num_buffer_ + 1);
- buff_box.setBig(boost_direction_, i_lab);
- BoxArray buff_ba(buff_box);
+ LabFrameDiags_[i]->buff_box_.setSmall(boost_direction_,
+ i_lab - num_buffer_ + 1);
+ LabFrameDiags_[i]->buff_box_.setBig(boost_direction_, i_lab);
+
+ BoxArray buff_ba(LabFrameDiags_[i]->buff_box_);
buff_ba.maxSize(max_box_size_);
DistributionMapping buff_dm(buff_ba);
- data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp_to_dump, 0) );
+ LabFrameDiags_[i]->data_buffer_.reset( new MultiFab(buff_ba,
+ buff_dm, ncomp_to_dump, 0) );
}
// ... reset particle buffer particles_buffer_[i]
if (WarpX::do_boosted_frame_particles)
- particles_buffer_[i].resize(mypc.nSpeciesBoostedFrameDiags());
+ LabFrameDiags_[i]->particles_buffer_.resize(
+ mypc.nSpeciesBoostedFrameDiags());
}
if (WarpX::do_boosted_frame_fields) {
const int ncomp = cell_centered_data->nComp();
const int start_comp = 0;
const bool interpolate = true;
- // Get slice in the boosted frame
- std::unique_ptr<MultiFab> slice = amrex::get_slice_data(boost_direction_,
- snapshots_[i].current_z_boost,
- *cell_centered_data, geom,
- start_comp, ncomp, interpolate);
-
- // transform it to the lab frame
- LorentzTransformZ(*slice, gamma_boost_, beta_boost_, ncomp);
- // Create a 2D box for the slice in the boosted frame
- Real dx = geom.CellSize(boost_direction_);
- int i_boost = (snapshots_[i].current_z_boost - geom.ProbLo(boost_direction_))/dx;
- Box slice_box = geom.Domain();
- slice_box.setSmall(boost_direction_, i_boost);
- slice_box.setBig(boost_direction_, i_boost);
- // Make it a BoxArray slice_ba
- BoxArray slice_ba(slice_box);
- slice_ba.maxSize(max_box_size_);
- // Create MultiFab tmp on slice_ba with data from slice
- MultiFab tmp(slice_ba, data_buffer_[i]->DistributionMap(), ncomp, 0);
- tmp.copy(*slice, 0, 0, ncomp);
-
- // Copy data from MultiFab tmp to MultiDab data_buffer[i]
- CopySlice(tmp, *data_buffer_[i], i_lab, map_actual_fields_to_dump);
+ // slice containing back-transformed data is generated only if t_lab != prev_t_lab and is re-used if multiple diags have the same z_lab,t_lab.
+ if (LabFrameDiags_[i]->t_lab != prev_t_lab ) {
+ if (slice)
+ {
+ slice.reset(new MultiFab);
+ slice.reset(nullptr);
+ }
+ slice = amrex::get_slice_data(boost_direction_,
+ LabFrameDiags_[i]->current_z_boost,
+ *cell_centered_data, geom,
+ start_comp, ncomp,
+ interpolate);
+ // Back-transform data to the lab-frame
+ LorentzTransformZ(*slice, gamma_boost_, beta_boost_, ncomp);
+ }
+ // Create a 2D box for the slice in the boosted frame
+ Real dx = geom.CellSize(boost_direction_);
+ int i_boost = ( LabFrameDiags_[i]->current_z_boost -
+ geom.ProbLo(boost_direction_))/dx;
+ //Box slice_box = geom.Domain();
+ Box slice_box = LabFrameDiags_[i]->buff_box_;
+ slice_box.setSmall(boost_direction_, i_boost);
+ slice_box.setBig(boost_direction_, i_boost);
+
+ // Make it a BoxArray slice_ba
+ BoxArray slice_ba(slice_box);
+ slice_ba.maxSize(max_box_size_);
+ tmp_slice_ptr = std::unique_ptr<MultiFab>(new MultiFab(slice_ba,
+ LabFrameDiags_[i]->data_buffer_->DistributionMap(),
+ ncomp, 0));
+
+ // slice is re-used if the t_lab of a diag is equal to
+ // that of the previous diag.
+ // Back-transformed data is copied from slice
+ // which has the dmap of the domain to
+ // tmp_slice_ptr which has the dmap of the
+ // data_buffer that stores the back-transformed data.
+ tmp_slice_ptr->copy(*slice, 0, 0, ncomp);
+ LabFrameDiags_[i]->AddDataToBuffer(*tmp_slice_ptr, i_lab,
+ map_actual_fields_to_dump);
+ tmp_slice_ptr.reset(new MultiFab);
+ tmp_slice_ptr.reset(nullptr);
}
if (WarpX::do_boosted_frame_particles) {
- mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_,
- old_z_boost, snapshots_[i].current_z_boost,
- t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]);
- }
+ if (LabFrameDiags_[i]->t_lab != prev_t_lab ) {
+ if (tmp_particle_buffer.size()>0)
+ {
+ tmp_particle_buffer.clear();
+ tmp_particle_buffer.shrink_to_fit();
+ }
+ tmp_particle_buffer.resize(mypc.nSpeciesBoostedFrameDiags());
+ mypc.GetLabFrameData( LabFrameDiags_[i]->file_name, i_lab,
+ boost_direction_, old_z_boost,
+ LabFrameDiags_[i]->current_z_boost,
+ t_boost, LabFrameDiags_[i]->t_lab, dt,
+ tmp_particle_buffer);
+ }
+ LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer,
+ mypc.nSpeciesBoostedFrameDiags());
+ }
- ++buff_counter_[i];
+ ++LabFrameDiags_[i]->buff_counter_;
+ prev_t_lab = LabFrameDiags_[i]->t_lab;
// If buffer full, write to disk.
- if (buff_counter_[i] == num_buffer_) {
+ if ( LabFrameDiags_[i]->buff_counter_ == num_buffer_) {
if (WarpX::do_boosted_frame_fields) {
#ifdef WARPX_USE_HDF5
- for (int comp = 0; comp < data_buffer_[i]->nComp(); ++comp)
- output_write_field(snapshots_[i].file_name, mesh_field_names[comp],
- *data_buffer_[i], comp);
+
+ Box buff_box = LabFrameDiags_[i]->buff_box_;
+ for (int comp = 0; comp < LabFrameDiags_[i]->data_buffer_->nComp(); ++comp)
+ output_write_field( LabFrameDiags_[i]->file_name,
+ mesh_field_names[comp],
+ *LabFrameDiags_[i]->data_buffer_, comp,
+ lbound(buff_box).x, lbound(buff_box).y,
+ lbound(buff_box).z);
#else
std::stringstream mesh_ss;
- mesh_ss << snapshots_[i].file_name << "/Level_0/" << Concatenate("buffer", i_lab, 5);
- VisMF::Write(*data_buffer_[i], mesh_ss.str());
+ mesh_ss << LabFrameDiags_[i]->file_name << "/Level_0/" <<
+ Concatenate("buffer", i_lab, 5);
+ VisMF::Write( (*LabFrameDiags_[i]->data_buffer_), mesh_ss.str());
#endif
}
@@ -788,24 +920,27 @@ writeLabFrameData(const MultiFab* cell_centered_data,
// Loop over species to be dumped to BFD
for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) {
// Get species name
- const std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(j)];
+ const std::string species_name = species_names[
+ mypc.mapSpeciesBoostedFrameDiags(j)];
#ifdef WARPX_USE_HDF5
// Write data to disk (HDF5)
- writeParticleDataHDF5(particles_buffer_[i][j],
- snapshots_[i].file_name,
+ writeParticleDataHDF5(LabFrameDiags_[i]->particles_buffer_[j],
+ LabFrameDiags_[i]->file_name,
species_name);
#else
std::stringstream part_ss;
- part_ss << snapshots_[i].file_name + "/" + species_name + "/";
+ part_ss << LabFrameDiags_[i]->file_name + "/" +
+ species_name + "/";
// Write data to disk (custom)
- writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab);
+ writeParticleData(LabFrameDiags_[i]->particles_buffer_[j],
+ part_ss.str(), i_lab);
#endif
}
- particles_buffer_[i].clear();
+ LabFrameDiags_[i]->particles_buffer_.clear();
}
- buff_counter_[i] = 0;
+ LabFrameDiags_[i]->buff_counter_ = 0;
}
}
@@ -823,7 +958,8 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat
Vector<long> particle_counts(ParallelDescriptor::NProcs(), 0);
Vector<long> particle_offsets(ParallelDescriptor::NProcs(), 0);
- ParallelAllGather::AllGather(np, particle_counts.data(), ParallelContext::CommunicatorAll());
+ ParallelAllGather::AllGather(np, particle_counts.data(),
+ ParallelContext::CommunicatorAll());
long total_np = 0;
for (int i = 0; i < ParallelDescriptor::NProcs(); ++i) {
@@ -854,7 +990,8 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat
output_write_particle_field(name, field_path,
pdata.GetRealData(k).data(),
particle_counts[ParallelDescriptor::MyProc()],
- particle_offsets[ParallelDescriptor::MyProc()] + old_np);
+ particle_offsets[ParallelDescriptor::MyProc()]
+ + old_np);
}
}
#endif
@@ -871,7 +1008,6 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata,
const int MyProc = ParallelDescriptor::MyProc();
auto np = pdata.GetRealData(DiagIdx::w).size();
-
if (np == 0) return;
field_name = name + Concatenate("w_", i_lab, 5) + "_" + std::to_string(MyProc);
@@ -917,14 +1053,14 @@ writeMetaData ()
BL_PROFILE("BoostedFrameDiagnostic::writeMetaData");
if (ParallelDescriptor::IOProcessor()) {
-
- if (!UtilCreateDirectory(WarpX::lab_data_directory, 0755))
- CreateDirectoryFailed(WarpX::lab_data_directory);
+ const std::string fullpath = WarpX::lab_data_directory + "/snapshots";
+ if (!UtilCreateDirectory(fullpath, 0755))
+ CreateDirectoryFailed(fullpath);
VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
std::ofstream HeaderFile;
HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
- std::string HeaderFileName(WarpX::lab_data_directory + "/Header");
+ std::string HeaderFileName(WarpX::lab_data_directory + "/snapshots/Header");
HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out |
std::ofstream::trunc |
std::ofstream::binary);
@@ -937,31 +1073,81 @@ writeMetaData ()
HeaderFile << dt_snapshots_lab_ << "\n";
HeaderFile << gamma_boost_ << "\n";
HeaderFile << beta_boost_ << "\n";
+
+ if (N_slice_snapshots_ > 0) {
+ const std::string fullpath_slice = WarpX::lab_data_directory + "/slices";
+ if (!UtilCreateDirectory(fullpath_slice, 0755))
+ CreateDirectoryFailed(fullpath_slice);
+
+ VisMF::IO_Buffer io_buffer_slice(VisMF::IO_Buffer_Size);
+ std::ofstream HeaderFile_slice;
+ HeaderFile_slice.rdbuf()->pubsetbuf(io_buffer_slice.dataPtr(),
+ io_buffer_slice.size());
+ std::string HeaderFileName_slice(WarpX::lab_data_directory+
+ "/slices/Header");
+ HeaderFile_slice.open(HeaderFileName_slice.c_str(),
+ std::ofstream::out |
+ std::ofstream::trunc |
+ std::ofstream::binary);
+
+ if (!HeaderFile_slice.good())
+ FileOpenFailed(HeaderFileName_slice);
+
+ HeaderFile_slice.precision(17);
+
+ HeaderFile_slice << N_slice_snapshots_ << "\n";
+ HeaderFile_slice << dt_slice_snapshots_lab_ << "\n";
+ HeaderFile_slice << gamma_boost_ << "\n";
+ HeaderFile_slice << beta_boost_ << "\n";
+
+ }
+
}
+
+
}
-BoostedFrameDiagnostic::LabSnapShot::
-LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab,
- IntVect prob_ncells_lab,
- int ncomp_to_dump,
- std::vector<std::string> mesh_field_names,
- int file_num_in,
- const BoostedFrameDiagnostic& bfd)
- : t_lab(t_lab_in),
- prob_domain_lab_(prob_domain_lab),
- prob_ncells_lab_(prob_ncells_lab),
- ncomp_to_dump_(ncomp_to_dump),
- mesh_field_names_(mesh_field_names),
- file_num(file_num_in),
- my_bfd(bfd)
+LabFrameSnapShot::
+LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in,
+ Real inv_beta_boost_in, Real dz_lab_in, RealBox prob_domain_lab,
+ IntVect prob_ncells_lab, int ncomp_to_dump,
+ std::vector<std::string> mesh_field_names,
+ amrex::RealBox diag_domain_lab, Box diag_box, int file_num_in)
{
- Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1);
- current_z_lab = 0.0;
- current_z_boost = 0.0;
- updateCurrentZPositions(t_boost, my_bfd.inv_gamma_boost_, my_bfd.inv_beta_boost_);
- initial_i = (current_z_lab - zmin_lab) / my_bfd.dz_lab_;
- file_name = Concatenate(WarpX::lab_data_directory + "/snapshot", file_num, 5);
+ t_lab = t_lab_in;
+ dz_lab_ = dz_lab_in;
+ inv_gamma_boost_ = inv_gamma_boost_in;
+ inv_beta_boost_ = inv_beta_boost_in;
+ prob_domain_lab_ = prob_domain_lab;
+ prob_ncells_lab_ = prob_ncells_lab;
+ diag_domain_lab_ = diag_domain_lab;
+ buff_box_ = diag_box;
+ ncomp_to_dump_ = ncomp_to_dump;
+ mesh_field_names_ = mesh_field_names;
+ file_num = file_num_in;
+ current_z_lab = 0.0;
+ current_z_boost = 0.0;
+ updateCurrentZPositions(t_boost, inv_gamma_boost_, inv_beta_boost_);
+ Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1);
+ initial_i = (current_z_lab - zmin_lab) / dz_lab_ ;
+ file_name = Concatenate(WarpX::lab_data_directory + "/snapshots/snapshot",
+ file_num, 5);
+ createLabFrameDirectories();
+ buff_counter_ = 0;
+ if (WarpX::do_boosted_frame_fields) data_buffer_.reset(nullptr);
+}
+
+void
+LabFrameDiag::
+updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta)
+{
+ current_z_boost = (t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta;
+ current_z_lab = (t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta;
+}
+void
+LabFrameDiag::
+createLabFrameDirectories() {
#ifdef WARPX_USE_HDF5
if (ParallelDescriptor::IOProcessor())
{
@@ -974,7 +1160,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab,
{
if (WarpX::do_boosted_frame_fields)
{
- for (int comp = 0; comp < ncomp_to_dump; ++comp) {
+ const auto lo = lbound(buff_box_);
+ for (int comp = 0; comp < ncomp_to_dump_; ++comp) {
output_create_field(file_name, mesh_field_names_[comp],
prob_ncells_lab_[0],
#if ( AMREX_SPACEDIM == 3 )
@@ -1036,19 +1223,12 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab,
#endif
ParallelDescriptor::Barrier();
- writeSnapShotHeader();
-}
-
-void
-BoostedFrameDiagnostic::LabSnapShot::
-updateCurrentZPositions(Real t_boost, Real inv_gamma, Real inv_beta) {
- current_z_boost = (t_lab*inv_gamma - t_boost)*PhysConst::c*inv_beta;
- current_z_lab = (t_lab - t_boost*inv_gamma)*PhysConst::c*inv_beta;
+ writeLabFrameHeader();
}
void
-BoostedFrameDiagnostic::LabSnapShot::
-writeSnapShotHeader() {
+LabFrameDiag::
+writeLabFrameHeader() {
#ifndef WARPX_USE_HDF5
if (ParallelDescriptor::IOProcessor()) {
VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
@@ -1072,17 +1252,17 @@ writeSnapShotHeader() {
<< prob_ncells_lab_[AMREX_SPACEDIM-1] <<'\n';
// Write domain physical boundaries
// domain lower bound
- HeaderFile << prob_domain_lab_.lo(0) << ' '
+ HeaderFile << diag_domain_lab_.lo(0) << ' '
#if ( AMREX_SPACEDIM==3 )
- << prob_domain_lab_.lo(1) << ' '
+ << diag_domain_lab_.lo(1) << ' '
#endif
- << prob_domain_lab_.lo(AMREX_SPACEDIM-1) <<'\n';
+ << diag_domain_lab_.lo(AMREX_SPACEDIM-1) <<'\n';
// domain higher bound
- HeaderFile << prob_domain_lab_.hi(0) << ' '
+ HeaderFile << diag_domain_lab_.hi(0) << ' '
#if ( AMREX_SPACEDIM==3 )
- << prob_domain_lab_.hi(1) << ' '
+ << diag_domain_lab_.hi(1) << ' '
#endif
- << prob_domain_lab_.hi(AMREX_SPACEDIM-1) <<'\n';
+ << diag_domain_lab_.hi(AMREX_SPACEDIM-1) <<'\n';
// List of fields dumped to file
for (int i=0; i<ncomp_to_dump_; i++)
{
@@ -1091,4 +1271,197 @@ writeSnapShotHeader() {
HeaderFile << "\n";
}
#endif
+
+}
+
+
+LabFrameSlice::
+LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in,
+ Real inv_beta_boost_in, Real dz_lab_in, RealBox prob_domain_lab,
+ IntVect prob_ncells_lab, int ncomp_to_dump,
+ std::vector<std::string> mesh_field_names,
+ RealBox diag_domain_lab, Box diag_box, int file_num_in,
+ amrex::Real cell_dx, amrex::Real cell_dy)
+{
+ t_lab = t_lab_in;
+ prob_domain_lab_ = prob_domain_lab;
+ prob_ncells_lab_ = prob_ncells_lab;
+ diag_domain_lab_ = diag_domain_lab;
+ buff_box_ = diag_box;
+ ncomp_to_dump_ = ncomp_to_dump;
+ mesh_field_names_ = mesh_field_names;
+ file_num = file_num_in;
+ current_z_lab = 0.0;
+ current_z_boost = 0.0;
+ updateCurrentZPositions(t_boost, inv_gamma_boost_, inv_beta_boost_);
+ Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1);
+ initial_i = (current_z_lab - zmin_lab)/dz_lab_;
+ file_name = Concatenate(WarpX::lab_data_directory+"/slices/slice",file_num,5);
+ createLabFrameDirectories();
+ buff_counter_ = 0;
+ dx_ = cell_dx;
+ dy_ = cell_dy;
+
+ if (WarpX::do_boosted_frame_fields) data_buffer_.reset(nullptr);
+}
+
+void
+LabFrameSnapShot::
+AddDataToBuffer( MultiFab& tmp, int k_lab,
+ amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump)
+{
+ const int ncomp_to_dump = map_actual_fields_to_dump.size();
+ MultiFab& buf = *data_buffer_;
+ for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+ Array4<Real> tmp_arr = tmp[mfi].array();
+ Array4<Real> buf_arr = buf[mfi].array();
+ // For 3D runs, tmp is a 2D (x,y) multifab that contains only
+ // slice to write to file
+ const Box& bx = mfi.tilebox();
+ const auto field_map_ptr = map_actual_fields_to_dump.dataPtr();
+ ParallelFor(bx, ncomp_to_dump,
+ [=] AMREX_GPU_DEVICE(int i, int j, int k, int n)
+ {
+ const int icomp = field_map_ptr[n];
+#if (AMREX_SPACEDIM == 3)
+ buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp);
+#else
+ buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp);
+#endif
+ }
+ );
+ }
+}
+
+
+void
+LabFrameSlice::
+AddDataToBuffer( MultiFab& tmp, int k_lab,
+ amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump)
+{
+ const int ncomp_to_dump = map_actual_fields_to_dump.size();
+ MultiFab& buf = *data_buffer_;
+ for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi)
+ {
+ Box& bx = buff_box_;
+ const Box& bx_bf = mfi.tilebox();
+ bx.setSmall(AMREX_SPACEDIM-1,bx_bf.smallEnd(AMREX_SPACEDIM-1));
+ bx.setBig(AMREX_SPACEDIM-1,bx_bf.bigEnd(AMREX_SPACEDIM-1));
+ Array4<Real> tmp_arr = tmp[mfi].array();
+ Array4<Real> buf_arr = buf[mfi].array();
+ const auto field_map_ptr = map_actual_fields_to_dump.dataPtr();
+ ParallelFor(bx, ncomp_to_dump,
+ [=] AMREX_GPU_DEVICE(int i, int j, int k, int n)
+ {
+ const int icomp = field_map_ptr[n];
+#if (AMREX_SPACEDIM == 3)
+ buf_arr(i,j,k_lab,n) = tmp_arr(i,j,k,icomp);
+#else
+ buf_arr(i,k_lab,k,n) = tmp_arr(i,j,k,icomp);
+#endif
+ });
+ }
+
+}
+
+
+void
+LabFrameSnapShot::
+AddPartDataToParticleBuffer(
+ Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer,
+ int nspeciesBoostedFrame) {
+ for (int isp = 0; isp < nspeciesBoostedFrame; ++isp) {
+ auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size();
+ if (np == 0) return;
+
+ particles_buffer_[isp].GetRealData(DiagIdx::w).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::w).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::w).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::w).end());
+
+ particles_buffer_[isp].GetRealData(DiagIdx::x).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::x).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::x).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::x).end());
+
+ particles_buffer_[isp].GetRealData(DiagIdx::y).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::y).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::y).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::y).end());
+
+ particles_buffer_[isp].GetRealData(DiagIdx::z).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::z).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::z).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::z).end());
+
+ particles_buffer_[isp].GetRealData(DiagIdx::ux).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::ux).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::ux).end());
+
+ particles_buffer_[isp].GetRealData(DiagIdx::uy).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::uy).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::uy).end());
+
+ particles_buffer_[isp].GetRealData(DiagIdx::uz).insert(
+ particles_buffer_[isp].GetRealData(DiagIdx::uz).end(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).begin(),
+ tmp_particle_buffer[isp].GetRealData(DiagIdx::uz).end());
+ }
+}
+
+void
+LabFrameSlice::
+AddPartDataToParticleBuffer(
+ Vector<WarpXParticleContainer::DiagnosticParticleData> tmp_particle_buffer,
+ int nSpeciesBoostedFrame) {
+ for (int isp = 0; isp < nSpeciesBoostedFrame; ++isp) {
+ auto np = tmp_particle_buffer[isp].GetRealData(DiagIdx::w).size();
+
+ if (np == 0) return;
+
+ auto const& wpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::w);
+ auto const& xpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::x);
+ auto const& ypc = tmp_particle_buffer[isp].GetRealData(DiagIdx::y);
+ auto const& zpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::z);
+ auto const& uxpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::ux);
+ auto const& uypc = tmp_particle_buffer[isp].GetRealData(DiagIdx::uy);
+ auto const& uzpc = tmp_particle_buffer[isp].GetRealData(DiagIdx::uz);
+
+ particles_buffer_[isp].resize(np);
+ auto wpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::w);
+ auto xpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::x);
+ auto ypc_buff = particles_buffer_[isp].GetRealData(DiagIdx::y);
+ auto zpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::z);
+ auto uxpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::ux);
+ auto uypc_buff = particles_buffer_[isp].GetRealData(DiagIdx::uy);
+ auto uzpc_buff = particles_buffer_[isp].GetRealData(DiagIdx::uz);
+
+
+ int partcounter = 0;
+ for (int i = 0; i < np; ++i)
+ {
+ if( xpc[i] >= (diag_domain_lab_.lo(0)-dx_) &&
+ xpc[i] <= (diag_domain_lab_.hi(0)+dx_) ) {
+ #if (AMREX_SPACEDIM == 3)
+ if( ypc[i] >= (diag_domain_lab_.lo(1)-dy_) &&
+ ypc[i] <= (diag_domain_lab_.hi(1) + dy_))
+ #endif
+ {
+ wpc_buff[partcounter] = wpc[i];
+ xpc_buff[partcounter] = xpc[i];
+ ypc_buff[partcounter] = ypc[i];
+ zpc_buff[partcounter] = zpc[i];
+ uxpc_buff[partcounter] = uxpc[i];
+ uypc_buff[partcounter] = uypc[i];
+ uzpc_buff[partcounter] = uzpc[i];
+ ++partcounter;
+ }
+ }
+ }
+
+ particles_buffer_[isp].resize(partcounter);
+
+ }
}
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 385993f78..78eaebfc5 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -86,17 +86,18 @@ WarpX::InitDiagnostics () {
const Real* current_lo = geom[0].ProbLo();
const Real* current_hi = geom[0].ProbHi();
Real dt_boost = dt[0];
-
// Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0
Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost );
Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost );
-
myBFD.reset(new BoostedFrameDiagnostic(zmin_lab,
zmax_lab,
moving_window_v, dt_snapshots_lab,
- num_snapshots_lab, gamma_boost,
- t_new[0], dt_boost,
- moving_window_dir, geom[0]));
+ num_snapshots_lab,
+ dt_slice_snapshots_lab,
+ num_slice_snapshots_lab,
+ gamma_boost, t_new[0], dt_boost,
+ moving_window_dir, geom[0],
+ slice_realbox));
}
}
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index c3ea3b763..d12a4dbff 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -18,6 +18,7 @@
#include <UpdateMomentumBoris.H>
#include <UpdateMomentumVay.H>
#include <UpdateMomentumBorisWithRadiationReaction.H>
+#include <UpdateMomentumHigueraCary.H>
using namespace amrex;
@@ -1011,12 +1012,12 @@ PhysicalParticleContainer::FieldGather (int lev,
const FArrayBox& byfab = By[pti];
const FArrayBox& bzfab = Bz[pti];
- Exp.assign(np,0.0);
- Eyp.assign(np,0.0);
- Ezp.assign(np,0.0);
- Bxp.assign(np,0.0);
- Byp.assign(np,0.0);
- Bzp.assign(np,0.0);
+ Exp.assign(np,WarpX::E_external[0]);
+ Eyp.assign(np,WarpX::E_external[1]);
+ Ezp.assign(np,WarpX::E_external[2]);
+ Bxp.assign(np,WarpX::B_external[0]);
+ Byp.assign(np,WarpX::B_external[1]);
+ Bzp.assign(np,WarpX::B_external[2]);
//
// copy data from particle container to temp arrays
@@ -1147,9 +1148,10 @@ PhysicalParticleContainer::Evolve (int lev,
exfab, eyfab, ezfab, bxfab, byfab, bzfab);
}
- Exp.assign(np,0.0);
- Eyp.assign(np,0.0);
- Ezp.assign(np,0.0);
+ Exp.assign(np,WarpX::E_external[0]);
+ Eyp.assign(np,WarpX::E_external[1]);
+ Ezp.assign(np,WarpX::E_external[2]);
+
Bxp.assign(np,WarpX::B_external[0]);
Byp.assign(np,WarpX::B_external[1]);
Bzp.assign(np,WarpX::B_external[2]);
@@ -1638,6 +1640,19 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
ux[i], uy[i], uz[i], dt );
}
);
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumHigueraCary( ux[i], uy[i], uz[i],
+ Ex[i], Ey[i], Ez[i], Bx[i],
+ By[i], Bz[i], qp, m, dt);
+ UpdatePosition( x[i], y[i], z[i],
+ ux[i], uy[i], uz[i], dt );
+ }
+ );
} else {
amrex::Abort("Unknown particle pusher");
};
@@ -1686,9 +1701,10 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
const FArrayBox& byfab = By[pti];
const FArrayBox& bzfab = Bz[pti];
- Exp.assign(np,0.0);
- Eyp.assign(np,0.0);
- Ezp.assign(np,0.0);
+ Exp.assign(np,WarpX::E_external[0]);
+ Eyp.assign(np,WarpX::E_external[1]);
+ Ezp.assign(np,WarpX::E_external[2]);
+
Bxp.assign(np,WarpX::B_external[0]);
Byp.assign(np,WarpX::B_external[1]);
Bzp.assign(np,WarpX::B_external[2]);
@@ -1764,6 +1780,13 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
qp, m, dt);
}
);
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
+ amrex::ParallelFor( pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ UpdateMomentumHigueraCary( ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ }
+ );
} else {
amrex::Abort("Unknown particle pusher");
};
diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package
index 45886702e..b033bfb57 100644
--- a/Source/Particles/Pusher/Make.package
+++ b/Source/Particles/Pusher/Make.package
@@ -2,6 +2,7 @@ CEXE_headers += GetAndSetPosition.H
CEXE_headers += UpdatePosition.H
CEXE_headers += UpdateMomentumBoris.H
CEXE_headers += UpdateMomentumVay.H
+CEXE_headers += UpdateMomentumHigueraCary.H
CEXE_headers += UpdatePositionPhoton.H
INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher
VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher
diff --git a/Source/Particles/Pusher/UpdateMomentumHigueraCary.H b/Source/Particles/Pusher/UpdateMomentumHigueraCary.H
new file mode 100644
index 000000000..51d7fd620
--- /dev/null
+++ b/Source/Particles/Pusher/UpdateMomentumHigueraCary.H
@@ -0,0 +1,59 @@
+#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_HIGUERACARY_H_
+#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_HIGUERACARY_H_
+
+#include "WarpXConst.H"
+#include <AMReX_FArrayBox.H>
+#include <AMReX_REAL.H>
+
+/**
+ * \brief Push the particle's positions over one timestep,
+ * given the value of its momenta `ux`, `uy`, `uz`
+ */
+
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void UpdateMomentumHigueraCary(
+ amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
+ const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
+ const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
+ const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+{
+ // Constants
+ const amrex::Real qmt = 0.5*q*dt/m;
+ constexpr amrex::Real invclight = 1./PhysConst::c;
+ constexpr amrex::Real invclightsq = 1./(PhysConst::c*PhysConst::c);
+ // Compute u_minus
+ const amrex::Real umx = ux + qmt*Ex;
+ const amrex::Real umy = uy + qmt*Ey;
+ const amrex::Real umz = uz + qmt*Ez;
+ // Compute gamma squared of u_minus
+ amrex::Real gamma = 1. + (umx*umx + umy*umy + umz*umz)*invclightsq;
+ // Compute beta and betam squared
+ const amrex::Real betax = qmt*Bx;
+ const amrex::Real betay = qmt*By;
+ const amrex::Real betaz = qmt*Bz;
+ const amrex::Real betam = betax*betax + betay*betay + betaz*betaz;
+ // Compute sigma
+ const amrex::Real sigma = gamma - betam;
+ // Get u*
+ const amrex::Real ust = (umx*betax + umy*betay + umz*betaz)*invclight;
+ // Get new gamma inversed
+ gamma = 1./std::sqrt(0.5*(sigma + std::sqrt(sigma*sigma + 4.*(betam + ust*ust)) ));
+ // Compute t
+ const amrex::Real tx = gamma*betax;
+ const amrex::Real ty = gamma*betay;
+ const amrex::Real tz = gamma*betaz;
+ // Compute s
+ const amrex::Real s = 1./(1.+(tx*tx + ty*ty + tz*tz));
+ // Compute um dot t
+ const amrex::Real umt = umx*tx + umy*ty + umz*tz;
+ // Compute u_plus
+ const amrex::Real upx = s*( umx + umt*tx + umy*tz - umz*ty );
+ const amrex::Real upy = s*( umy + umt*ty + umz*tx - umx*tz );
+ const amrex::Real upz = s*( umz + umt*tz + umx*ty - umy*tx );
+ // Get new u
+ ux = upx + qmt*Ex + upy*tz - upz*ty;
+ uy = upy + qmt*Ey + upz*tx - upx*tz;
+ uz = upz + qmt*Ez + upx*ty - upy*tx;
+}
+
+#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_HIGUERACARY_H_
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index 24d4ac24c..535ffec6f 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -14,6 +14,7 @@
#include <UpdateMomentumBoris.H>
#include <UpdateMomentumVay.H>
#include <UpdateMomentumBorisWithRadiationReaction.H>
+#include <UpdateMomentumHigueraCary.H>
using namespace amrex;
@@ -391,9 +392,9 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
const FArrayBox& byfab = By[pti];
const FArrayBox& bzfab = Bz[pti];
- Exp.assign(np,0.0);
- Eyp.assign(np,0.0);
- Ezp.assign(np,0.0);
+ Exp.assign(np,WarpX::E_external[0]);
+ Eyp.assign(np,WarpX::E_external[1]);
+ Ezp.assign(np,WarpX::E_external[2]);
Bxp.assign(np,WarpX::B_external[0]);
Byp.assign(np,WarpX::B_external[1]);
Bzp.assign(np,WarpX::B_external[2]);
@@ -463,6 +464,13 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
q, m, dt);
}
);
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
+ amrex::ParallelFor( pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ UpdateMomentumHigueraCary( uxpp[i], uypp[i], uzpp[i],
+ Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ }
+ );
} else {
amrex::Abort("Unknown particle pusher");
};
diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H
index 54c721abf..269171a5f 100644
--- a/Source/Utils/WarpXAlgorithmSelection.H
+++ b/Source/Utils/WarpXAlgorithmSelection.H
@@ -14,7 +14,8 @@ struct MaxwellSolverAlgo {
struct ParticlePusherAlgo {
enum {
Boris = 0,
- Vay = 1
+ Vay = 1,
+ HigueraCary = 2
};
};
diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp
index 8a6ff6dbf..be9cdd118 100644
--- a/Source/Utils/WarpXAlgorithmSelection.cpp
+++ b/Source/Utils/WarpXAlgorithmSelection.cpp
@@ -18,6 +18,7 @@ const std::map<std::string, int> maxwell_solver_algo_to_int = {
const std::map<std::string, int> particle_pusher_algo_to_int = {
{"boris", ParticlePusherAlgo::Boris },
{"vay", ParticlePusherAlgo::Vay },
+ {"higuera", ParticlePusherAlgo::HigueraCary },
{"default", ParticlePusherAlgo::Boris }
};
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 0da1cf350..29b89686e 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -75,6 +75,7 @@ public:
// External fields
static amrex::Vector<amrex::Real> B_external;
+ static amrex::Vector<amrex::Real> E_external;
// Algorithms
static long current_deposition_algo;
@@ -269,6 +270,9 @@ public:
void SliceGenerationForDiagnostics ();
void WriteSlicePlotFile () const;
void ClearSliceMultiFabs ();
+ static int num_slice_snapshots_lab;
+ static amrex::Real dt_slice_snapshots_lab;
+ amrex::RealBox getSliceRealBox() const {return slice_realbox;}
// these should be private, but can't due to Cuda limitations
static void ComputeDivB (amrex::MultiFab& divB, int dcomp,
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 5a51d7d13..d94541f17 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -25,6 +25,7 @@
using namespace amrex;
Vector<Real> WarpX::B_external(3, 0.0);
+Vector<Real> WarpX::E_external(3, 0.0);
int WarpX::do_moving_window = 0;
int WarpX::moving_window_dir = -1;
@@ -67,6 +68,9 @@ Real WarpX::dt_snapshots_lab = std::numeric_limits<Real>::lowest();
bool WarpX::do_boosted_frame_fields = true;
bool WarpX::do_boosted_frame_particles = true;
+int WarpX::num_slice_snapshots_lab = 0;
+Real WarpX::dt_slice_snapshots_lab;
+
bool WarpX::do_dynamic_scheduling = true;
int WarpX::do_subcycling = 0;
@@ -76,9 +80,9 @@ IntVect WarpX::Bx_nodal_flag(1,0,0);
IntVect WarpX::By_nodal_flag(0,1,0);
IntVect WarpX::Bz_nodal_flag(0,0,1);
#elif (AMREX_SPACEDIM == 2)
-IntVect WarpX::Bx_nodal_flag(1,0); // x is the first dimension to AMReX
-IntVect WarpX::By_nodal_flag(0,0); // y is the missing dimension to 2D AMReX
-IntVect WarpX::Bz_nodal_flag(0,1); // z is the second dimension to 2D AMReX
+IntVect WarpX::Bx_nodal_flag(1,0);// x is the first dimension to AMReX
+IntVect WarpX::By_nodal_flag(0,0);// y is the missing dimension to 2D AMReX
+IntVect WarpX::Bz_nodal_flag(0,1);// z is the second dimension to 2D AMReX
#endif
#if (AMREX_SPACEDIM == 3)
@@ -86,9 +90,9 @@ IntVect WarpX::Ex_nodal_flag(0,1,1);
IntVect WarpX::Ey_nodal_flag(1,0,1);
IntVect WarpX::Ez_nodal_flag(1,1,0);
#elif (AMREX_SPACEDIM == 2)
-IntVect WarpX::Ex_nodal_flag(0,1); // x is the first dimension to AMReX
-IntVect WarpX::Ey_nodal_flag(1,1); // y is the missing dimension to 2D AMReX
-IntVect WarpX::Ez_nodal_flag(1,0); // z is the second dimension to 2D AMReX
+IntVect WarpX::Ex_nodal_flag(0,1);// x is the first dimension to AMReX
+IntVect WarpX::Ey_nodal_flag(1,1);// y is the missing dimension to 2D AMReX
+IntVect WarpX::Ez_nodal_flag(1,0);// z is the second dimension to 2D AMReX
#endif
#if (AMREX_SPACEDIM == 3)
@@ -96,9 +100,9 @@ IntVect WarpX::jx_nodal_flag(0,1,1);
IntVect WarpX::jy_nodal_flag(1,0,1);
IntVect WarpX::jz_nodal_flag(1,1,0);
#elif (AMREX_SPACEDIM == 2)
-IntVect WarpX::jx_nodal_flag(0,1); // x is the first dimension to AMReX
-IntVect WarpX::jy_nodal_flag(1,1); // y is the missing dimension to 2D AMReX
-IntVect WarpX::jz_nodal_flag(1,0); // z is the second dimension to 2D AMReX
+IntVect WarpX::jx_nodal_flag(0,1);// x is the first dimension to AMReX
+IntVect WarpX::jy_nodal_flag(1,1);// y is the missing dimension to 2D AMReX
+IntVect WarpX::jz_nodal_flag(1,0);// z is the second dimension to 2D AMReX
#endif
IntVect WarpX::filter_npass_each_dir(1);
@@ -253,13 +257,13 @@ void
WarpX::ReadParameters ()
{
{
- ParmParse pp; // Traditionally, max_step and stop_time do not have prefix.
+ ParmParse pp;// Traditionally, max_step and stop_time do not have prefix.
pp.query("max_step", max_step);
pp.query("stop_time", stop_time);
}
{
- ParmParse pp("amr"); // Traditionally, these have prefix, amr.
+ ParmParse pp("amr");// Traditionally, these have prefix, amr.
pp.query("check_file", check_file);
pp.query("check_int", check_int);
@@ -291,6 +295,7 @@ WarpX::ReadParameters ()
zmax_plasma_to_compute_max_step);
pp.queryarr("B_external", B_external);
+ pp.queryarr("E_external", E_external);
pp.query("do_moving_window", do_moving_window);
if (do_moving_window)
@@ -591,6 +596,15 @@ WarpX::ReadParameters ()
}
}
+ if (do_boosted_frame_diagnostic) {
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0,
+ "gamma_boost must be > 1 to use the boost frame diagnostic");
+ pp.query("num_slice_snapshots_lab", num_slice_snapshots_lab);
+ if (num_slice_snapshots_lab > 0) {
+ pp.get("dt_slice_snapshots_lab", dt_slice_snapshots_lab );
+ }
+ }
+
}
}
@@ -915,8 +929,8 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
RealVect cdx_vect(cdx[0], cdx[2]);
#endif
// Get the cell-centered box, with guard cells
- BoxArray realspace_ba = cba; // Copy box
- realspace_ba.enclosedCells().grow(ngE); // cell-centered + guard cells
+ BoxArray realspace_ba = cba;// Copy box
+ realspace_ba.enclosedCells().grow(ngE);// cell-centered + guard cells
// Define spectral solver
spectral_solver_cp[lev].reset( new SpectralSolver( realspace_ba, dm,
nox_fft, noy_fft, noz_fft, do_nodal, cdx_vect, dt[lev] ) );
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-12 20:28:27 +0000