diff options
21 files changed, 1226 insertions, 302 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 019b907b5..e77ca553e 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -392,6 +392,10 @@ Particle initialization It only works if `<species>.do_qed = 1`. Enables non-linear Breit-Wheeler process for this species. **Implementation of this feature is in progress. It requires to compile with QED=TRUE** +* ``warpx.E_external`` & ``warpx.B_external`` (list of `float`) optional (default `0.0`) + Two separate parameters which add a uniform E-field or B-field to each particle + which is then added to the field values gathered from the grid in the + PIC cycle. Laser initialization -------------------- @@ -623,6 +627,7 @@ Numerics and algorithms - ``boris``: Boris pusher. - ``vay``: Vay pusher (see `Vay, Phys. Plasmas (2008) <https://aip.scitation.org/doi/10.1063/1.2837054>`__) + - ``higuera``: Higuera-Cary pusher (see `Higuera and Cary, Phys. Plasmas (2017) <https://aip.scitation.org/doi/10.1063/1.4979989>`__) If ``algo.particle_pusher`` is not specified, ``boris`` is the default. @@ -823,16 +828,37 @@ Diagnostics and output ``warpx.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz part_per_cell``. * ``slice.dom_lo`` and ``slice.dom_hi`` (`2 floats in 2D`, `3 floats in 3D`; in meters similar to the units of the simulation box.) - The extent of the slice are defined by the co-ordinates of the lower corner (``slice.dom_lo``) and upper corner (``slice.dom_hi``). The slice could be 1D, 2D, or 3D, aligned with the co-ordinate axes and the first axis of the coordinates is x. For example: if for a 3D simulation, an x-z slice is to be extracted at y = 0.0, then the y-value of slice.dom_lo and slice.dom_hi must be equal to 0.0 + The extent of the slice are defined by the co-ordinates of the lower + corner (``slice.dom_lo``) and upper corner (``slice.dom_hi``). + The slice could be 1D, 2D, or 3D, aligned with the co-ordinate axes + and the first axis of the coordinates is x. For example: if for a + 3D simulation, an x-z slice is to be extracted at y = 0.0, + then the y-value of slice.dom_lo and slice.dom_hi must be equal to 0.0 * ``slice.coarsening_ratio`` (`2 integers in 2D`, `3 integers in 3D`; default `1`) The coarsening ratio input must be greater than 0. Default is 1 in all directions. - In the directions that is reduced, i.e., for an x-z slice in 3D, the reduced y-dimension has a default coarsening ratio equal to 1. + In the directions that is reduced, i.e., for an x-z slice in 3D, + the reduced y-dimension has a default coarsening ratio equal to 1. * ``slice.plot_int`` (`integer`) The number of PIC cycles inbetween two consecutive data dumps for the slice. Use a negative number to disable slice generation and slice data dumping. +* ``slice.num_slice_snapshots_lab`` (`integer`) + Only used when ``warpx.do_boosted_frame_diagnostic`` is ``1``. + The number of back-transformed field and particle data that + will be written for the reduced domain defined by ``slice.dom_lo`` + and ``slice.dom_hi``. Note that the 'slice' is a reduced + diagnostic which could be 1D, 2D, or 3D, aligned with the co-ordinate axes. + These slices can be visualized using read_raw_data.py and the HDF5 format can + be visualized using the h5py library. Please see the documentation on visualization + for further details. + +* ``slice.dt_slice_snapshots_lab`` (`float`, in seconds) + Only used when ``warpx.do_boosted_frame_diagnostic`` is ``1``. + The time interval between the back-transformed reduced diagnostics (where this + time interval is expressed in the laboratory frame). + Checkpoints and restart ----------------------- WarpX supports checkpoints/restart via AMReX. diff --git a/Docs/source/visualization/yt.rst b/Docs/source/visualization/yt.rst index 2e0a13fcb..e2f8b377f 100644 --- a/Docs/source/visualization/yt.rst +++ b/Docs/source/visualization/yt.rst @@ -124,6 +124,50 @@ Alternatively, the main commands can be found in our example jupyter notebook for postprocessing :download:`Visualization.ipynb<../../../Tools/Visualization.ipynb>`. +The full back-transformed diagnostics of the entire domain is written in lab_frame_data/snapshots/ and the back-transformed diagnostics of the reduced domain is written to lab_frame_data/slices/ +For instance : To plot the ``Ez`` field along the z-direction at the center of the 3D-domain of the full back-transformed diagnostics for the entire 3D domain : + +:: + iteration = 0 + field = 'Ez' + snapshot = './lab_frame_data/snapshots/' + 'snapshot' + str(iteration).zfill(5) + header = './lab_frame_data/snapshots/Header' + allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) # Read field data + F = allrd[field] + plt.plot(F[F.shape[0]//2,F.shape[1]//2-1,:]) + +Similarly, the back-transformed diagnostics on a reduced domain (1D line, 2D slice, 3D reduced diagnostic) can also be visualized using read_raw_data.py. For instance -- let us say that the user-input is an x-z slice at y=y_mid of the domain, then, to plot ``Ez`` of the x-z slice along the z-direction at the center of the slice : + +:: + iteration = 0 + field = 'Ez' + snapshot = './lab_frame_data/slices/' + 'slice' + str(iteration).zfill(5) + header = './lab_frame_data/slices/Header' + allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) # Read field data + F_RD = allrd[field] + plt.plot(F_RD[F_RD.shape[0]//2,0,:]) + + +Note that, in the above snippet, we compare the 0th cell of the reduced diagnostic with F.shape[1]//2-1. For an x-z slice at y=y-mid of the domain, two cells are extracted at the center to ensure that the data format is HDF5 compliant. Let us consider that the domain consists of four cells in the y-dimension : [0,1,2,3], Then the 2D slice would contain the data that corresponds to [1,2]. That is the 0th cell of the reduced diagnostic corresponds to ny/2-1, (where, ny is the number of cells in the y-dimension). + +If the back-transformed diagnostics are written in the HDF5 format (This can be done by compiling WarpX with USE_HDF5=TRUE), then the full domain snapshot and reduced domain diagnostics can be visualized using h5py : + +:: + import matplotlib.pyplot as plt + import h5py + f1 = h5py.File('lab_frame_data/snapshots/snapshot00000', 'r') + nx1 = f1['Ez'].shape[0] + ny1 = f1['Ez'].shape[1] + nz1 = f1['Ez'].shape[2] + plt.plot(f1['Ez'][nx1//2,ny1//2-1,:]) + + f2 = h5py.File('lab_frame_data/slices/slice00000', 'r') + nx2 = f2['Ez'].shape[0] + ny2 = f2['Ez'].shape[1] + nz2 = f2['Ez'].shape[2] + plt.figure() + plt.plot(f2['Ez'][nx2//2,0,:]) + Further information ------------------- diff --git a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py index 3dcf40ac8..f28272897 100755 --- a/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py +++ b/Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py @@ -19,8 +19,8 @@ import read_raw_data yt.funcs.mylog.setLevel(0) # Read data from back-transformed diagnostics -snapshot = './lab_frame_data/snapshot00001' -header = './lab_frame_data/Header' +snapshot = './lab_frame_data/snapshots/snapshot00001' +header = './lab_frame_data/snapshots/Header' allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) z = np.mean( read_raw_data.get_particle_field(snapshot, 'beam', 'z') ) w = np.std ( read_raw_data.get_particle_field(snapshot, 'beam', 'x') ) diff --git a/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py b/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py new file mode 100755 index 000000000..3b8e7aa76 --- /dev/null +++ b/Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py @@ -0,0 +1,36 @@ +#! /usr/bin/env python + +''' +Analysis script of a WarpX simulation in a boosted frame. + +The simulation runs in a boosted frame, and the analysis is done in the lab +frame, i.e., on the back-transformed diagnostics for the full 3D simulation and +an x-z slice at y=y_center. The field-data, Ez, along z, at (x_center,y_center,:) is compared +between the full back-transformed diagnostic and the reduced diagnostic (i.e., x-z slice) . +''' + +import numpy as np +import read_raw_data + +# Read data from back-transformed diagnostics of entire domain +snapshot = './lab_frame_data/snapshots/snapshot00002' +header = './lab_frame_data/snapshots/Header' +allrd, info = read_raw_data.read_lab_snapshot(snapshot, header) +F = allrd['Ez'] +print("F.shape ", F.shape) +F_1D = np.squeeze(F[F.shape[0]//2,F.shape[1]//2,:]) + + +# Read data from reduced back-transformed diagnostics (i.e. slice) +snapshot_slice = './lab_frame_data/slices/slice00002' +header_slice = './lab_frame_data/slices/Header' +allrd, info = read_raw_data.read_lab_snapshot(snapshot_slice, header_slice) +Fs = allrd['Ez'] +print("Fs.shape", Fs.shape) +Fs_1D = np.squeeze(Fs[Fs.shape[0]//2,1,:]) + +error = np.max(np.abs(Fs_1D - F_1D)) / np.max(np.abs(F_1D)) + +# Print error and assert small error +print("relative error: " + str(error)) +assert( error < 1E-15 ) diff --git a/Examples/Modules/boosted_diags/inputs.3d.slice b/Examples/Modules/boosted_diags/inputs.3d.slice new file mode 100644 index 000000000..fdd150949 --- /dev/null +++ b/Examples/Modules/boosted_diags/inputs.3d.slice @@ -0,0 +1,111 @@ +warpx.zmax_plasma_to_compute_max_step = 0.0031 + +amr.n_cell = 32 32 64 +amr.max_grid_size = 64 +amr.blocking_factor = 32 +amr.max_level = 0 +amr.plot_int = 10000 + +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 1 1 0 # Is periodic? +geometry.prob_lo = -128.e-6 -128.e-6 -40.e-6 +geometry.prob_hi = 128.e-6 128.e-6 0.96e-6 + +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = standard +algo.particle_pusher = vay +algo.maxwell_fdtd_solver = ckc +interpolation.nox = 3 +interpolation.noy = 3 +interpolation.noz = 3 +warpx.use_filter = 1 +warpx.cfl = 1. +warpx.do_pml = 0 + +warpx.do_moving_window = 1 +warpx.moving_window_dir = z +warpx.moving_window_v = 1.0 # in units of the speed of light +warpx.serialize_ics = 1 + +warpx.gamma_boost = 10. +warpx.boost_direction = z +warpx.do_boosted_frame_diagnostic = 1 +warpx.num_snapshots_lab = 4 +warpx.dz_snapshots_lab = 0.001 +warpx.boosted_frame_diag_fields= Ex Ey Ez By rho + +particles.nspecies = 3 +particles.species_names = electrons ions beam +particles.use_fdtd_nci_corr = 1 + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = NUniformPerCell +electrons.num_particles_per_cell_each_dim = 1 1 1 +electrons.momentum_distribution_type = "gaussian" +electrons.xmin = -120.e-6 +electrons.xmax = 120.e-6 +electrons.ymin = -120.e-6 +electrons.ymax = 120.e-6 +electrons.zmin = 0. +electrons.zmax = .003 +electrons.profile = constant +electrons.density = 3.5e24 +electrons.do_continuous_injection = 1 +electrons.do_boosted_frame_diags = 1 + +ions.charge = q_e +ions.mass = m_p +ions.injection_style = NUniformPerCell +ions.num_particles_per_cell_each_dim = 1 1 1 +ions.momentum_distribution_type = "gaussian" +ions.xmin = -120.e-6 +ions.xmax = 120.e-6 +ions.ymin = -120.e-6 +ions.ymax = 120.e-6 +ions.zmin = 0. +ions.zmax = .003 +ions.profile = constant +ions.density = 3.5e24 +ions.do_continuous_injection = 1 +ions.do_boosted_frame_diags = 1 + +beam.charge = -q_e +beam.mass = m_e +beam.injection_style = "gaussian_beam" +beam.x_rms = 1.e-6 +beam.y_rms = 1.e-6 +beam.z_rms = .2e-6 +beam.x_m = 0. +beam.y_m = 0. +beam.z_m = -20.e-6 +beam.npart = 1000 +beam.q_tot = -1.e-14 +beam.momentum_distribution_type = "gaussian" +beam.ux_m = 0.0 +beam.uy_m = 0.0 +beam.uz_m = 200000. +beam.ux_th = .2 +beam.uy_th = .2 +beam.uz_th = 20. + +lasers.nlasers = 1 +lasers.names = laser1 +laser1.profile = Gaussian +laser1.position = 0. 0. -0.1e-6 # This point is on the laser plane +laser1.direction = 0. 0. 1. # The plane normal direction +laser1.polarization = 0. 1. 0. # The main polarization vector +laser1.e_max = 2.e12 # Maximum amplitude of the laser field (in V/m) +laser1.profile_waist = 45.e-6 # The waist of the laser (in meters) +laser1.profile_duration = 20.e-15 # The duration of the laser (in seconds) +laser1.profile_t_peak = 40.e-15 # The time at which the laser reaches its peak (in seconds) +laser1.profile_focal_distance = 0.5e-3 # Focal distance from the antenna (in meters) +laser1.wavelength = 0.81e-6 # The wavelength of the laser (in meters) + +slice.dom_lo = -128.e-6 0.0 -40.e-6 +slice.dom_hi = 128.e-6 0.0 0.96e-6 +slice.coarsening_ratio = 1 1 1 +slice.plot_int = -1 +slice.num_slice_snapshots_lab = 4 +slice.dt_slice_snapshots_lab = 3.3356409519815207e-12 diff --git a/Examples/Tests/Langmuir/inputs.rt b/Examples/Tests/Langmuir/inputs.rt index be5a730c2..28602f758 100644 --- a/Examples/Tests/Langmuir/inputs.rt +++ b/Examples/Tests/Langmuir/inputs.rt @@ -26,6 +26,7 @@ warpx.verbose = 1 # Algorithms algo.field_gathering = standard +algo.particle_pusher = "higuera" # Interpolation interpolation.nox = 1 diff --git a/Examples/Tests/particle_pusher/analysis_pusher.py b/Examples/Tests/particle_pusher/analysis_pusher.py new file mode 100755 index 000000000..0d9fc24c5 --- /dev/null +++ b/Examples/Tests/particle_pusher/analysis_pusher.py @@ -0,0 +1,30 @@ +#! /usr/bin/env python + +# This script tests the particle pusher (HC) +# using a force-free field, +# in which position x should remain 0. +# An initial velocity Vy corresponding to +# Lorentz factor = 20 is used. +# Bz is fixed at 1 T. +# Ex = -Vy*Bz. + +# Possible errors: +# Boris: 2321.3958529 +# Vay: 0.00010467 +# HC: 0.00011403 +# tolerance: 0.001 +# Possible running time: ~ 4.0 s + +import sys +import yt + +tolerance = 0.001 + +filename = sys.argv[1] +ds = yt.load( filename ) +ad = ds.all_data() +x = ad['particle_position_x'].to_ndarray() + +print('error = ', abs(x)) +print('tolerance = ', tolerance) +assert(abs(x) < tolerance) diff --git a/Examples/Tests/particle_pusher/inputs b/Examples/Tests/particle_pusher/inputs new file mode 100755 index 000000000..45c075f51 --- /dev/null +++ b/Examples/Tests/particle_pusher/inputs @@ -0,0 +1,44 @@ +# Maximum number of time steps +max_step = 10000 + +# number of grid points +amr.n_cell = 8 8 8 + +# Maximum level in hierarchy (disable mesh refinement) +amr.max_level = 0 + +# How often to write plotfiles. "<= 0" means no plotfiles. +amr.plot_int = 10000 + +# Geometry +geometry.coord_sys = 0 # Cartesian +geometry.is_periodic = 1 1 1 # yes +geometry.prob_lo = -2.077023075927835e+07 -2.077023075927835e+07 -2.077023075927835e+07 +geometry.prob_hi = 2.077023075927835e+07 2.077023075927835e+07 2.077023075927835e+07 + +# PML +warpx.do_pml = 0 + +# Algorithms +algo.charge_deposition = standard +algo.field_gathering = standard +algo.particle_pusher = "higuera" + +# CFL +warpx.cfl = 1.0 + +# particles +particles.nspecies = 1 +particles.species_names = positron +positron.charge = 1.0 +positron.mass = 1.0 +positron.injection_style = "SingleParticle" +positron.single_particle_pos = 0.0 0.0 0.0 +positron.single_particle_vel = 0.0 19.974984355438178 0.0 +positron.single_particle_weight = 0.0 +warpx.plot_raw_fields = 0 + +# External fields +# Ex is set to be Ex = -vy*Bz +warpx.B_external = 0.0 0.0 1.0 +warpx.E_external = -2.994174829214179e+08 0.0 0.0 diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 401c777b6..a3521cf6a 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -134,6 +134,23 @@ doComparison = 0 aux1File = Tools/read_raw_data.py analysisRoutine = Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py +[Boost_3Dbacktransformed_CheckReducedDiagnostics] +buildDir = . +inputFile = Examples/Modules/boosted_diags/inputs.3d.slice +dim = 3 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +compareParticles = 0 +doComparison = 0 +aux1File = Tools/read_raw_data.py +analysisRoutine = Examples/Modules/boosted_diags/analysis_3Dbacktransformed_diag.py + [nci_corrector] buildDir = . inputFile = Examples/Modules/nci_corrector/inputs2d @@ -677,3 +694,17 @@ compileTest = 0 doVis = 0 compareParticles = 0 analysisRoutine = Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py + +[particle_pusher] +buildDir = . +inputFile = Examples/Tests/particle_pusher/inputs +dim = 3 +restartTest = 0 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Tests/particle_pusher/analysis_pusher.py diff --git a/Regression/prepare_file_travis.py b/Regression/prepare_file_travis.py index 088b97b04..6c3be105f 100644 --- a/Regression/prepare_file_travis.py +++ b/Regression/prepare_file_travis.py @@ -41,15 +41,15 @@ text = re.sub( 'numprocs = \d+', 'numprocs = 1', text) text = re.sub( 'numthreads = \d+', 'numthreads = 1', text) # Remove Python test (does not compile) -text = re.sub( '\[Python_Langmuir\]\n(.+\n)*\n', '', text) +text = re.sub( '\[Python_Langmuir\]\n(.+\n)*', '', text) # Remove Langmuir_x/y/z test (too long; not that useful) -text = re.sub( '\[Langmuir_[xyz]\]\n(.+\n)*\n', '', text) +text = re.sub( '\[Langmuir_[xyz]\]\n(.+\n)*', '', text) # Remove tests that do not have the right dimension if dim is not None: print('Selecting tests with dim = %s' %dim) - text = re.sub('\[.+\n(.+\n)*dim = [^%s]\n(.+\n)*\n' %dim, '', text) + text = re.sub('\[.+\n(.+\n)*dim = [^%s]\n(.+\n)*' %dim, '', text) # Remove or keep QED tests according to 'qed' variable if qed is not None: 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] ) ); |