diff options
Diffstat (limited to 'Source')
57 files changed, 961 insertions, 961 deletions
diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index 3a09259b0..6c8eb15db 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -14,24 +14,24 @@ #include "WarpXConst.H" /// -/// BoostedFrameDiagnostic is for handling IO when running in a boosted +/// 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 +/// 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 + /// 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 + /// 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_; @@ -42,7 +42,7 @@ class BoostedFrameDiagnostic { int ncomp_to_dump_; std::vector<std::string> mesh_field_names_; - + int file_num; int initial_i; const BoostedFrameDiagnostic& my_bfd; @@ -54,12 +54,12 @@ class BoostedFrameDiagnostic { 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. + /// 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); @@ -87,17 +87,17 @@ class BoostedFrameDiagnostic { amrex::Vector<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 + // - i is the back-transformed snapshot number // - j is the species number amrex::Vector<amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> > particles_buffer_; int num_buffer_ = 256; int max_box_size_ = 256; - // buff_counter_[i] is the number of z slices in data_buffer_[i] + // buff_counter_[i] is the number of z slices in data_buffer_[i] // for snapshot number i. amrex::Vector<int> buff_counter_; amrex::Vector<LabSnapShot> snapshots_; - + void writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const int i_lab); @@ -106,44 +106,44 @@ class BoostedFrameDiagnostic { const std::string& name, const std::string& species_name); #endif public: - - BoostedFrameDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab, + + 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); - + void Flush(const amrex::Geometry& geom); - + void writeLabFrameData(const amrex::MultiFab* cell_centered_data, const MultiParticleContainer& mypc, const amrex::Geometry& geom, const amrex::Real t_boost, const amrex::Real dt); - + void writeMetaData(); private: // Map field names and component number in cell_centered_data std::map<std::string, int> possible_fields_to_dump = { - {"Ex" , 0}, - {"Ey" , 1}, - {"Ez" , 2}, - {"Bx" , 3}, - {"By" , 4}, - {"Bz" , 5}, - {"jx" , 6}, - {"jy" , 7}, - {"jz" , 8}, + {"Ex" , 0}, + {"Ey" , 1}, + {"Ez" , 2}, + {"Bx" , 3}, + {"By" , 4}, + {"Bz" , 5}, + {"jx" , 6}, + {"jy" , 7}, + {"jz" , 8}, {"rho", 9} }; - // maps field index in data_buffer_[i] -> cell_centered_data for + // maps field index in data_buffer_[i] -> cell_centered_data for // snapshots i. By default, all fields in cell_centered_data are dumped. // Needs to be amrex::Vector because used in a ParallelFor kernel. amrex::Gpu::ManagedDeviceVector<int> map_actual_fields_to_dump; // Name of fields to dump. By default, all fields in cell_centered_data. // Needed for file headers only. - std::vector<std::string> mesh_field_names = {"Ex", "Ey", "Ez", - "Bx", "By", "Bz", + std::vector<std::string> mesh_field_names = {"Ex", "Ey", "Ez", + "Bx", "By", "Bz", "jx", "jy", "jz", "rho"}; int ncomp_to_dump = 10; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index ad4f50806..f1d437509 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -17,7 +17,7 @@ using namespace amrex; namespace { 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. Should be run only by the root process. @@ -32,7 +32,7 @@ namespace } /* - Writes a single string attribute to the given group. + Writes a single string attribute to the given group. Should only be called by the root process. */ void write_string_attribute(hid_t& group, const std::string& key, const std::string& val) @@ -52,7 +52,7 @@ namespace } /* - Writes a single double attribute to the given group. + Writes a single double attribute to the given group. Should only be called by the root process. */ void write_double_attribute(hid_t& group, const std::string& key, const double val) @@ -109,7 +109,7 @@ namespace */ void output_create_field(const std::string& file_path, const std::string& field_path, const unsigned nx, const unsigned ny, const unsigned nz) - { + { BL_PROFILE("output_create_field"); // Open the output. @@ -121,7 +121,7 @@ namespace hsize_t dims[3] = {nx, nz}; #endif hid_t grid_space = H5Screate_simple(AMREX_SPACEDIM, dims, NULL); - + // Create the dataset. hid_t dataset = H5Dcreate(file, field_path.c_str(), H5T_IEEE_F64LE, grid_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); @@ -141,7 +141,7 @@ namespace /* Creates a group associated with a single particle species. Should be run by all processes collectively. - */ + */ void output_create_species_group(const std::string& file_path, const std::string& species_name) { MPI_Comm comm = MPI_COMM_WORLD; @@ -152,10 +152,10 @@ namespace // Create the file access prop list. hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS); H5Pset_fapl_mpio(pa_plist, comm, info); - + // Open the output. hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist); - + hid_t group = H5Gcreate(file, species_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Gclose(group); @@ -169,7 +169,7 @@ namespace */ long output_resize_particle_field(const std::string& file_path, const std::string& field_path, const long num_to_add) - { + { BL_PROFILE("output_resize_particle_field"); // Open the output. @@ -187,7 +187,7 @@ namespace hsize_t new_size[1]; new_size[0] = dims[0] + num_to_add; status = H5Dset_extent (dataset, new_size); - + if (status < 0) { amrex::Abort("Error: set extent filed on dataset " @@ -208,7 +208,7 @@ namespace */ void output_write_particle_field(const std::string& file_path, const std::string& field_path, const Real* data_ptr, const long count, const long index) - { + { BL_PROFILE("output_write_particle_field"); MPI_Comm comm = MPI_COMM_WORLD; @@ -219,15 +219,15 @@ namespace // Create the file access prop list. hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS); H5Pset_fapl_mpio(pa_plist, comm, info); - + // Open the output. hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist); - + int RANK = 1; hsize_t offset[1]; hsize_t dims[1]; herr_t status; - + hid_t dataset = H5Dopen (file, field_path.c_str(), H5P_DEFAULT); // Make sure the dataset is there. @@ -236,9 +236,9 @@ namespace amrex::Abort("Error on rank " + std::to_string(mpi_rank) + ". Count not find dataset " + field_path + "\n"); } - + hid_t filespace = H5Dget_space (dataset); - + offset[0] = index; dims[0] = count; @@ -247,23 +247,23 @@ namespace H5Pset_dxpl_mpio(collective_plist, H5FD_MPIO_INDEPENDENT); if (count > 0) { - + /* Define memory space */ hid_t memspace = H5Screate_simple (RANK, dims, NULL); - + status = H5Sselect_hyperslab (filespace, H5S_SELECT_SET, offset, NULL, dims, NULL); - + if (status < 0) { amrex::Abort("Error on rank " + std::to_string(ParallelDescriptor::MyProc()) + " could not select hyperslab.\n"); } - + /* Write the data to the extended portion of dataset */ status = H5Dwrite(dataset, H5T_NATIVE_DOUBLE, memspace, filespace, collective_plist, data_ptr); - + if (status < 0) { amrex::Abort("Error on rank " + std::to_string(ParallelDescriptor::MyProc()) + @@ -272,9 +272,9 @@ namespace status = H5Sclose (memspace); } - + ParallelDescriptor::Barrier(); - + // Close resources. H5Pclose(collective_plist); H5Sclose(filespace); @@ -282,13 +282,13 @@ namespace H5Fclose(file); H5Pclose(pa_plist); } - + /* Creates an extendible dataset, suitable for storing particle data. Should be run on all ranks collectively. */ void output_create_particle_field(const std::string& file_path, const std::string& field_path) - { + { BL_PROFILE("output_create_particle_field"); MPI_Comm comm = MPI_COMM_WORLD; @@ -299,7 +299,7 @@ namespace // Create the file access prop list. hid_t pa_plist = H5Pcreate(H5P_FILE_ACCESS); H5Pset_fapl_mpio(pa_plist, comm, info); - + // Open the output. hid_t file = H5Fopen(file_path.c_str(), H5F_ACC_RDWR, pa_plist); @@ -307,7 +307,7 @@ namespace hsize_t dims[1] = {0}; hsize_t maxdims[1] = {H5S_UNLIMITED}; hsize_t chunk_dims[2] = {4}; - + hid_t dataspace = H5Screate_simple (RANK, dims, maxdims); // Enable chunking @@ -316,7 +316,7 @@ namespace hid_t dataset = H5Dcreate2 (file, field_path.c_str(), H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT, prop, H5P_DEFAULT); - + if (dataset < 0) { amrex::Abort("Error: could not create dataset. H5 returned " @@ -329,7 +329,7 @@ namespace H5Sclose(dataspace); H5Fclose(file); } - + /* Write the only component in the multifab to the dataset given by field_name. Uses hdf5-parallel. @@ -361,7 +361,7 @@ namespace amrex::Abort("Error on rank " + std::to_string(mpi_rank) + ". Count not find dataset " + field_path + "\n"); } - + // Grab the dataspace of the field dataset from file. hid_t file_dataspace = H5Dget_space(dataset); @@ -388,9 +388,9 @@ namespace const Box& box = mfi.validbox(); const int *lo_vec = box.loVect(); const int *hi_vec = box.hiVect(); - + transposed_data.resize(box.numPts(), 0.0); - + // Set slab offset and shape. for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { @@ -399,7 +399,7 @@ namespace slab_offsets[idim] = lo_vec[idim]; slab_dims[idim] = hi_vec[idim] - lo_vec[idim] + 1; } - + int cnt = 0; AMREX_D_TERM( for (int i = lo_vec[0]; i <= hi_vec[0]; ++i), @@ -433,7 +433,7 @@ namespace } ParallelDescriptor::Barrier(); - + // Close HDF5 resources. H5Pclose(collective_plist); H5Sclose(file_dataspace); @@ -447,7 +447,7 @@ namespace namespace { void - CopySlice(MultiFab& tmp, MultiFab& buf, int k_lab, + 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(); @@ -458,11 +458,11 @@ namespace 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 + // 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(); + 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) { @@ -488,13 +488,13 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) for (MFIter mfi(data, TilingIfNotGPU()); mfi.isValid(); ++mfi) { const Box& tile_box = mfi.tilebox(); Array4< Real > arr = data[mfi].array(); - // arr(x,y,z,comp) where 0->9 comps are + // arr(x,y,z,comp) where 0->9 comps are // Ex Ey Ez Bx By Bz jx jy jz rho Real clight = PhysConst::c; ParallelFor(tile_box, [=] AMREX_GPU_DEVICE (int i, int j, int k) { - // Transform the transverse E and B fields. Note that ez and bz are not + // 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)); @@ -523,8 +523,8 @@ 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 gamma_boost, Real t_boost, Real dt_boost, + Real dt_snapshots_lab, int N_snapshots, + Real gamma_boost, Real t_boost, Real dt_boost, int boost_direction, const Geometry& geom) : gamma_boost_(gamma_boost), dt_snapshots_lab_(dt_snapshots_lab), @@ -537,11 +537,11 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, AMREX_ALWAYS_ASSERT(WarpX::do_boosted_frame_fields or WarpX::do_boosted_frame_particles); - + inv_gamma_boost_ = 1.0 / gamma_boost_; beta_boost_ = std::sqrt(1.0 - inv_gamma_boost_*inv_gamma_boost_); inv_beta_boost_ = 1.0 / beta_boost_; - + dz_lab_ = PhysConst::c * dt_boost_ * inv_beta_boost_ * inv_gamma_boost_; inv_dz_lab_ = 1.0 / dz_lab_; int Nz_lab = static_cast<unsigned>((zmax_lab - zmin_lab) * inv_dz_lab_); @@ -563,9 +563,9 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, std::vector<std::string> user_fields_to_dump; ParmParse pp("warpx"); bool do_user_fields; - do_user_fields = pp.queryarr("boosted_frame_diag_fields", + do_user_fields = pp.queryarr("boosted_frame_diag_fields", user_fields_to_dump); - // If user specifies fields to dump, overwrite ncomp_to_dump, + // If user specifies fields to dump, overwrite ncomp_to_dump, // map_actual_fields_to_dump and mesh_field_names. for (int i = 0; i < 10; ++i) map_actual_fields_to_dump.push_back(i); if (do_user_fields){ @@ -588,7 +588,7 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, 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, + LabSnapShot snapshot(t_lab, t_boost, prob_domain_lab, prob_ncells_lab, ncomp_to_dump, mesh_field_names, i, *this); snapshots_.push_back(snapshot); @@ -602,54 +602,54 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, void BoostedFrameDiagnostic::Flush(const Geometry& geom) { BL_PROFILE("BoostedFrameDiagnostic::Flush"); - + VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); - + // Loop over BFD snapshots for (int i = 0; i < N_snapshots_; ++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_; - + if (buff_counter_[i] != 0) { if (WarpX::do_boosted_frame_fields) { const BoxArray& ba = data_buffer_[i]->boxArray(); const int hi = ba[0].bigEnd(boost_direction_); const int lo = hi - buff_counter_[i] + 1; - + Box buff_box = geom.Domain(); buff_box.setSmall(boost_direction_, lo); buff_box.setBig(boost_direction_, hi); - + BoxArray buff_ba(buff_box); buff_ba.maxSize(max_box_size_); DistributionMapping buff_dm(buff_ba); - + const int ncomp = data_buffer_[i]->nComp(); - + MultiFab tmp(buff_ba, buff_dm, ncomp, 0); - + tmp.copy(*data_buffer_[i], 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); -#else +#else std::stringstream ss; ss << snapshots_[i].file_name << "/Level_0/" << Concatenate("buffer", i_lab, 5); VisMF::Write(tmp, ss.str()); #endif } - + if (WarpX::do_boosted_frame_particles) { // Loop over species to be dumped to BFD for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { // Get species name - std::string species_name = + std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; #ifdef WARPX_USE_HDF5 // Dump species data @@ -691,7 +691,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // Loop over snapshots for (int i = 0; i < N_snapshots_; ++i) { - + // Get updated z position of snapshot const Real old_z_boost = snapshots_[i].current_z_boost; snapshots_[i].updateCurrentZPositions(t_boost, @@ -700,15 +700,15 @@ writeLabFrameData(const MultiFab* cell_centered_data, Real zmin_lab = snapshots_[i].prob_domain_lab_.lo(AMREX_SPACEDIM-1); Real zmax_lab = snapshots_[i].prob_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; - // Get z index of data_buffer_ (i.e. in the lab frame) where - // simulation domain (t', [zmin',zmax']), back-transformed to lab + // 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_; @@ -725,7 +725,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp_to_dump, 0) ); } // ... reset particle buffer particles_buffer_[i] - if (WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nSpeciesBoostedFrameDiags()); } @@ -738,7 +738,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, 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 @@ -753,7 +753,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // 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); } @@ -766,7 +766,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, ++buff_counter_[i]; - + // If buffer full, write to disk. if (buff_counter_[i] == num_buffer_) { @@ -781,10 +781,10 @@ writeLabFrameData(const MultiFab* cell_centered_data, VisMF::Write(*data_buffer_[i], mesh_ss.str()); #endif } - + if (WarpX::do_boosted_frame_particles) { // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { // Get species name const std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; #ifdef WARPX_USE_HDF5 @@ -800,14 +800,14 @@ writeLabFrameData(const MultiFab* cell_centered_data, // Write data to disk (custom) writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif - } + } particles_buffer_[i].clear(); } buff_counter_[i] = 0; } } - - VisMF::SetHeaderVersion(current_version); + + VisMF::SetHeaderVersion(current_version); } #ifdef WARPX_USE_HDF5 @@ -817,12 +817,12 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat const std::string& name, const std::string& species_name) { auto np = pdata.GetRealData(DiagIdx::w).size(); - + Vector<long> particle_counts(ParallelDescriptor::NProcs(), 0); Vector<long> particle_offsets(ParallelDescriptor::NProcs(), 0); - + ParallelAllGather::AllGather(np, particle_counts.data(), ParallelContext::CommunicatorAll()); - + long total_np = 0; for (int i = 0; i < ParallelDescriptor::NProcs(); ++i) { particle_offsets[i] = total_np; @@ -830,7 +830,7 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat } if (total_np == 0) return; - + long old_np = 0; if (ParallelDescriptor::IOProcessor()) { @@ -844,7 +844,7 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat // Note, this has the effect of an MPI Barrier between the above resize operation // and the below write. ParallelDescriptor::ReduceLongMax(old_np); - + // Write data here for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k) { @@ -853,7 +853,7 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat pdata.GetRealData(k).data(), particle_counts[ParallelDescriptor::MyProc()], particle_offsets[ParallelDescriptor::MyProc()] + old_np); - } + } } #endif @@ -863,7 +863,7 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const int i_lab) { BL_PROFILE("BoostedFrameDiagnostic::writeParticleData"); - + std::string field_name; std::ofstream ofs; @@ -880,27 +880,27 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, field_name = name + Concatenate("x_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); writeRealData(pdata.GetRealData(DiagIdx::x).data(), np, ofs); - ofs.close(); + ofs.close(); field_name = name + Concatenate("y_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); writeRealData(pdata.GetRealData(DiagIdx::y).data(), np, ofs); - ofs.close(); + ofs.close(); field_name = name + Concatenate("z_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); writeRealData(pdata.GetRealData(DiagIdx::z).data(), np, ofs); - ofs.close(); - + ofs.close(); + field_name = name + Concatenate("ux_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); writeRealData(pdata.GetRealData(DiagIdx::ux).data(), np, ofs); - ofs.close(); + ofs.close(); field_name = name + Concatenate("uy_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); writeRealData(pdata.GetRealData(DiagIdx::uy).data(), np, ofs); - ofs.close(); + ofs.close(); field_name = name + Concatenate("uz_", i_lab, 5) + "_" + std::to_string(MyProc); ofs.open(field_name.c_str(), std::ios::out|std::ios::binary); @@ -910,12 +910,12 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, void BoostedFrameDiagnostic:: -writeMetaData () +writeMetaData () { BL_PROFILE("BoostedFrameDiagnostic::writeMetaData"); if (ParallelDescriptor::IOProcessor()) { - + if (!UtilCreateDirectory(WarpX::lab_data_directory, 0755)) CreateDirectoryFailed(WarpX::lab_data_directory); @@ -928,22 +928,22 @@ writeMetaData () std::ofstream::binary); if(!HeaderFile.good()) FileOpenFailed(HeaderFileName); - + HeaderFile.precision(17); - + HeaderFile << N_snapshots_ << "\n"; - HeaderFile << dt_snapshots_lab_ << "\n"; + HeaderFile << dt_snapshots_lab_ << "\n"; HeaderFile << gamma_boost_ << "\n"; HeaderFile << beta_boost_ << "\n"; } } BoostedFrameDiagnostic::LabSnapShot:: -LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, - IntVect prob_ncells_lab, +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, + int file_num_in, const BoostedFrameDiagnostic& bfd) : t_lab(t_lab_in), prob_domain_lab_(prob_domain_lab), @@ -967,7 +967,7 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, } ParallelDescriptor::Barrier(); - + if (ParallelDescriptor::IOProcessor()) { if (WarpX::do_boosted_frame_fields) @@ -986,7 +986,7 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, } ParallelDescriptor::Barrier(); - + if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); @@ -994,7 +994,7 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { // Loop over species to be dumped to BFD - std::string species_name = + std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k) @@ -1003,10 +1003,10 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, output_create_particle_field(file_name, field_path); } } - } -#else + } +#else if (ParallelDescriptor::IOProcessor()) { - + if (!UtilCreateDirectory(file_name, 0755)) CreateDirectoryFailed(file_name); @@ -1018,13 +1018,13 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, } auto & mypc = WarpX::GetInstance().GetPartContainer(); - const std::vector<std::string> species_names = mypc.GetSpeciesNames(); - + const std::vector<std::string> species_names = mypc.GetSpeciesNames(); + const std::string particles_prefix = "particle"; // Loop over species to be dumped to BFD for(int i = 0; i < mypc.nSpeciesBoostedFrameDiags(); ++i) { // Get species name - std::string species_name = + std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(i)]; const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) @@ -1058,9 +1058,9 @@ writeSnapShotHeader() { std::ofstream::binary); if(!HeaderFile.good()) FileOpenFailed(HeaderFileName); - + HeaderFile.precision(17); - + HeaderFile << t_lab << "\n"; // Write domain number of cells HeaderFile << prob_ncells_lab_[0] << ' ' diff --git a/Source/Diagnostics/ElectrostaticIO.cpp b/Source/Diagnostics/ElectrostaticIO.cpp index a8dffac2b..75be2a24c 100644 --- a/Source/Diagnostics/ElectrostaticIO.cpp +++ b/Source/Diagnostics/ElectrostaticIO.cpp @@ -114,7 +114,7 @@ WritePlotFileES (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, particle_varnames.push_back("Bz"); Vector<std::string> int_names; - + mypc->Checkpoint(plotfilename, particle_varnames, int_names); WriteJobInfo(plotfilename); diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 8bfa45a59..741d3f50d 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -22,8 +22,8 @@ RigidInjectedParticleContainer::ReadHeader (std::istream& is) int zinject_plane_tmp; is >> zinject_plane_tmp; zinject_plane_levels.push_back(zinject_plane_tmp); - WarpX::GotoNextLine(is); - } + WarpX::GotoNextLine(is); + } for (int i = 0; i < nlevs; ++i) { @@ -31,7 +31,7 @@ RigidInjectedParticleContainer::ReadHeader (std::istream& is) is >> done_injecting_tmp; done_injecting.push_back(done_injecting_tmp); WarpX::GotoNextLine(is); - } + } } void @@ -78,7 +78,7 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const { for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - auto& pc = allcontainers[i]; + auto& pc = allcontainers[i]; if (pc->plot_species) { Vector<std::string> real_names; @@ -90,19 +90,19 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); - + real_names.push_back("Ex"); real_names.push_back("Ey"); real_names.push_back("Ez"); - + real_names.push_back("Bx"); real_names.push_back("By"); real_names.push_back("Bz"); - + #ifdef WARPX_DIM_RZ real_names.push_back("theta"); #endif - + if(pc->do_field_ionization){ int_names.push_back("ionization_level"); // int_flags specifies, for each integer attribs, whether it is @@ -134,7 +134,7 @@ MultiParticleContainer::Restart (const std::string& dir) } void -MultiParticleContainer::ReadHeader (std::istream& is) +MultiParticleContainer::ReadHeader (std::istream& is) { for (auto& pc : allcontainers) { pc->ReadHeader(is); @@ -149,7 +149,7 @@ MultiParticleContainer::WriteHeader (std::ostream& os) const } } -// Particle momentum is defined as gamma*velocity, which is neither +// Particle momentum is defined as gamma*velocity, which is neither // SI mass*gamma*velocity nor normalized gamma*velocity/c. // This converts momentum to SI units (or vice-versa) to write SI data // to file. diff --git a/Source/Diagnostics/SliceDiagnostic.H b/Source/Diagnostics/SliceDiagnostic.H index 31eea83be..1b9ca3967 100644 --- a/Source/Diagnostics/SliceDiagnostic.H +++ b/Source/Diagnostics/SliceDiagnostic.H @@ -16,26 +16,26 @@ -std::unique_ptr<amrex::MultiFab> CreateSlice( const amrex::MultiFab& mf, - const amrex::Vector<amrex::Geometry> &dom_geom, - amrex::RealBox &slice_realbox, +std::unique_ptr<amrex::MultiFab> CreateSlice( const amrex::MultiFab& mf, + const amrex::Vector<amrex::Geometry> &dom_geom, + amrex::RealBox &slice_realbox, amrex::IntVect &slice_cr_ratio ); -void CheckSliceInput( const amrex::RealBox real_box, - amrex::RealBox &slice_cc_nd_box, amrex::RealBox &slice_realbox, +void CheckSliceInput( const amrex::RealBox real_box, + amrex::RealBox &slice_cc_nd_box, amrex::RealBox &slice_realbox, amrex::IntVect &slice_cr_ratio, amrex::Vector<amrex::Geometry> dom_geom, - amrex::IntVect const SliceType, amrex::IntVect &slice_lo, + amrex::IntVect const SliceType, amrex::IntVect &slice_lo, amrex::IntVect &slice_hi, amrex::IntVect &interp_lo); -void InterpolateSliceValues( amrex::MultiFab& smf, - amrex::IntVect interp_lo, amrex::RealBox slice_realbox, - amrex::Vector<amrex::Geometry> geom, int ncomp, int nghost, - amrex::IntVect slice_lo, amrex::IntVect slice_hi, +void InterpolateSliceValues( amrex::MultiFab& smf, + amrex::IntVect interp_lo, amrex::RealBox slice_realbox, + amrex::Vector<amrex::Geometry> geom, int ncomp, int nghost, + amrex::IntVect slice_lo, amrex::IntVect slice_hi, amrex::IntVect SliceType, const amrex::RealBox real_box); -void InterpolateLo( const amrex::Box& bx, amrex::FArrayBox &fabox, - amrex::IntVect slice_lo, amrex::Vector<amrex::Geometry> geom, - int idir, amrex::IntVect IndType, amrex::RealBox slice_realbox, +void InterpolateLo( const amrex::Box& bx, amrex::FArrayBox &fabox, + amrex::IntVect slice_lo, amrex::Vector<amrex::Geometry> geom, + int idir, amrex::IntVect IndType, amrex::RealBox slice_realbox, int srccomp, int ncomp, int nghost, const amrex::RealBox real_box); #endif diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp index 0ec528e32..79f03985b 100644 --- a/Source/Diagnostics/SliceDiagnostic.cpp +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -10,23 +10,23 @@ using namespace amrex; /* \brief * The functions creates the slice for diagnostics based on the user-input. - * The slice can be 1D, 2D, or 3D and it inherts the index type of the underlying data. - * The implementation assumes that the slice is aligned with the coordinate axes. - * The input parameters are modified if the user-input does not comply with requirements of coarsenability or if the slice extent is not contained within the simulation domain. - * First a slice multifab (smf) with cell size equal to that of the simulation grid is created such that it extends from slice.dim_lo to slice.dim_hi and shares the same index space as the source multifab (mf) + * The slice can be 1D, 2D, or 3D and it inherts the index type of the underlying data. + * The implementation assumes that the slice is aligned with the coordinate axes. + * The input parameters are modified if the user-input does not comply with requirements of coarsenability or if the slice extent is not contained within the simulation domain. + * First a slice multifab (smf) with cell size equal to that of the simulation grid is created such that it extends from slice.dim_lo to slice.dim_hi and shares the same index space as the source multifab (mf) * The values are copied from src mf to dst smf using amrex::ParallelCopy - * If interpolation is required, then on the smf, using data points stored in the ghost cells, the data in interpolated. - * If coarsening is required, then a coarse slice multifab is generated (cs_mf) and the - * values of the refined slice (smf) is averaged down to obtain the coarse slice. - * \param mf is the source multifab containing the field data - * \param dom_geom is the geometry of the domain and used in the function to obtain the - * CellSize of the underlying grid. - * \param slice_realbox defines the extent of the slice - * \param slice_cr_ratio provides the coarsening ratio for diagnostics + * If interpolation is required, then on the smf, using data points stored in the ghost cells, the data in interpolated. + * If coarsening is required, then a coarse slice multifab is generated (cs_mf) and the + * values of the refined slice (smf) is averaged down to obtain the coarse slice. + * \param mf is the source multifab containing the field data + * \param dom_geom is the geometry of the domain and used in the function to obtain the + * CellSize of the underlying grid. + * \param slice_realbox defines the extent of the slice + * \param slice_cr_ratio provides the coarsening ratio for diagnostics */ -std::unique_ptr<MultiFab> -CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, +std::unique_ptr<MultiFab> +CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, RealBox &slice_realbox, IntVect &slice_cr_ratio ) { std::unique_ptr<MultiFab> smf; @@ -36,10 +36,10 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, int nghost = 1; int nlevels = dom_geom.size(); int ncomp = (mf).nComp(); - - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nlevels==1, + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nlevels==1, "Slice diagnostics does not work with mesh refinement yet (TO DO)."); - + const auto conversionType = (mf).ixType(); IntVect SliceType(AMREX_D_DECL(0,0,0)); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim ) @@ -50,7 +50,7 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, const RealBox& real_box = dom_geom[0].ProbDomain(); RealBox slice_cc_nd_box; int slice_grid_size = 32; - + bool interpolate = false; bool coarsen = false; @@ -60,45 +60,45 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, IntVect interp_lo(AMREX_D_DECL(0,0,0)); CheckSliceInput(real_box, slice_cc_nd_box, slice_realbox, slice_cr_ratio, - dom_geom, SliceType, slice_lo, + dom_geom, SliceType, slice_lo, slice_hi, interp_lo); int configuration_dim = 0; // Determine if interpolation is required and number of cells in slice // - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + // Flag for interpolation if required // if ( interp_lo[idim] == 1) { - interpolate = 1; + interpolate = 1; } // For the case when a dimension is reduced // - if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { + if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { slice_ncells[idim] = 1; } - else { - slice_ncells[idim] = ( slice_hi[idim] - slice_lo[idim] + 1 ) + else { + slice_ncells[idim] = ( slice_hi[idim] - slice_lo[idim] + 1 ) / slice_cr_ratio[idim]; int refined_ncells = slice_hi[idim] - slice_lo[idim] + 1 ; - if ( slice_cr_ratio[idim] > 1) { + if ( slice_cr_ratio[idim] > 1) { coarsen = true; // modify slice_grid_size if >= refines_cells // if ( slice_grid_size >= refined_ncells ) { slice_grid_size = refined_ncells - 1; } - + } configuration_dim += 1; } } if (configuration_dim==1) { amrex::Warning("The slice configuration is 1D and cannot be visualized using yt."); - } + } // Slice generation with index type inheritance // Box slice(slice_lo, slice_hi); - + Vector<BoxArray> sba(1); sba[0].define(slice); sba[0].maxSize(slice_grid_size); @@ -106,7 +106,7 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, // Distribution mapping for slice can be different from that of domain // Vector<DistributionMapping> sdmap(1); sdmap[0] = DistributionMapping{sba[0]}; - + smf.reset(new MultiFab(amrex::convert(sba[0],SliceType), sdmap[0], ncomp, nghost)); @@ -115,12 +115,12 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, smf->ParallelCopy(mf, 0, 0, ncomp,nghost,nghost); // inteprolate if required on refined slice // - if (interpolate == 1 ) { - InterpolateSliceValues( *smf, interp_lo, slice_cc_nd_box, dom_geom, + if (interpolate == 1 ) { + InterpolateSliceValues( *smf, interp_lo, slice_cc_nd_box, dom_geom, ncomp, nghost, slice_lo, slice_hi, SliceType, real_box); } - + if (coarsen == false) { return smf; } @@ -131,14 +131,14 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, AMREX_ALWAYS_ASSERT(crse_ba[0].size() == sba[0].size()); - cs_mf.reset( new MultiFab(amrex::convert(crse_ba[0],SliceType), + cs_mf.reset( new MultiFab(amrex::convert(crse_ba[0],SliceType), sdmap[0], ncomp,nghost)); MultiFab& mfSrc = *smf; MultiFab& mfDst = *cs_mf; MFIter mfi_dst(mfDst); - for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { + for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { Array4<Real const> const& Src_fabox = mfSrc.const_array(mfi); @@ -150,7 +150,7 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, IntVect cctype(AMREX_D_DECL(0,0,0)); if( SliceType==cctype ) { - amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, + amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio); } IntVect ndtype(AMREX_D_DECL(1,1,1)); @@ -159,11 +159,11 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, scomp, ncomp, slice_cr_ratio); } if( SliceType == WarpX::Ex_nodal_flag ) { - amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 0); } if( SliceType == WarpX::Ey_nodal_flag) { - amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 1); } if( SliceType == WarpX::Ez_nodal_flag ) { @@ -171,19 +171,19 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, scomp, ncomp, slice_cr_ratio, 2); } if( SliceType == WarpX::Bx_nodal_flag) { - amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 0); } if( SliceType == WarpX::By_nodal_flag ) { - amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 1); } if( SliceType == WarpX::Bz_nodal_flag ) { - amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, ncomp, slice_cr_ratio, 2); } - if ( mfi_dst.isValid() ) { + if ( mfi_dst.isValid() ) { ++mfi_dst; } @@ -197,66 +197,66 @@ CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, /* \brief - * This function modifies the slice input parameters under certain conditions. - * The coarsening ratio, slice_cr_ratio is modified if the input is not an exponent of 2. - * for example, if the coarsening ratio is 3, 5 or 6, which is not an exponent of 2, - * then the value of coarsening ratio is modified to the nearest exponent of 2. - * The default value for coarsening ratio is 1. - * slice_realbox.lo and slice_realbox.hi are set equal to the simulation domain lo and hi - * if for the user-input for the slice lo and hi coordinates are outside the domain. - * If the slice_realbox.lo and slice_realbox.hi coordinates do not align with the data + * This function modifies the slice input parameters under certain conditions. + * The coarsening ratio, slice_cr_ratio is modified if the input is not an exponent of 2. + * for example, if the coarsening ratio is 3, 5 or 6, which is not an exponent of 2, + * then the value of coarsening ratio is modified to the nearest exponent of 2. + * The default value for coarsening ratio is 1. + * slice_realbox.lo and slice_realbox.hi are set equal to the simulation domain lo and hi + * if for the user-input for the slice lo and hi coordinates are outside the domain. + * If the slice_realbox.lo and slice_realbox.hi coordinates do not align with the data * points and the number of cells in that dimension is greater than 1, and if the extent of * the slice in that dimension is not coarsenable, then the value lo and hi coordinates are * shifted to the nearest coarsenable point to include some extra data points in the slice. * If slice_realbox.lo==slice_realbox.hi, then that dimension has only one cell and no - * modifications are made to the value. If the lo and hi do not align with a data point, - * then it is flagged for interpolation. - * \param real_box a Real box defined for the underlying domain. - * \param slice_realbox a Real box for defining the slice dimension. - * \param slice_cc_nd_box a Real box for defining the modified lo and hi of the slice + * modifications are made to the value. If the lo and hi do not align with a data point, + * then it is flagged for interpolation. + * \param real_box a Real box defined for the underlying domain. + * \param slice_realbox a Real box for defining the slice dimension. + * \param slice_cc_nd_box a Real box for defining the modified lo and hi of the slice * such that the coordinates align with the underlying data points. - * If the dimension is reduced to have only one cell, the slice_realbox is not modified and * instead the values are interpolated to the coordinate from the nearest data points. - * \param slice_cr_ratio contains values of the coarsening ratio which may be modified - * if the input values do not satisfy coarsenability conditions. + * If the dimension is reduced to have only one cell, the slice_realbox is not modified and * instead the values are interpolated to the coordinate from the nearest data points. + * \param slice_cr_ratio contains values of the coarsening ratio which may be modified + * if the input values do not satisfy coarsenability conditions. * \param slice_lo and slice_hi are the index values of the slice - * \param interp_lo are set to 0 or 1 if they are flagged for interpolation. - * The slice shares the same index space as that of the simulation domain. + * \param interp_lo are set to 0 or 1 if they are flagged for interpolation. + * The slice shares the same index space as that of the simulation domain. */ -void -CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, - RealBox &slice_realbox, IntVect &slice_cr_ratio, - Vector<Geometry> dom_geom, IntVect const SliceType, +void +CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, + RealBox &slice_realbox, IntVect &slice_cr_ratio, + Vector<Geometry> dom_geom, IntVect const SliceType, IntVect &slice_lo, IntVect &slice_hi, IntVect &interp_lo) { - + IntVect slice_lo2(AMREX_D_DECL(0,0,0)); - for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) - { + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { // Modify coarsening ratio if the input value is not an exponent of 2 for AMR // - if ( slice_cr_ratio[idim] > 0 ) { + if ( slice_cr_ratio[idim] > 0 ) { int log_cr_ratio = floor ( log2( double(slice_cr_ratio[idim]))); slice_cr_ratio[idim] = exp2( double(log_cr_ratio) ); } - + //// Default coarsening ratio is 1 // // Modify lo if input is out of bounds // - if ( slice_realbox.lo(idim) < real_box.lo(idim) ) { + if ( slice_realbox.lo(idim) < real_box.lo(idim) ) { slice_realbox.setLo( idim, real_box.lo(idim)); - amrex::Print() << " slice lo is out of bounds. " << - " Modified it in dimension " << idim << + amrex::Print() << " slice lo is out of bounds. " << + " Modified it in dimension " << idim << " to be aligned with the domain box\n"; - } - + } + // Modify hi if input in out od bounds // - if ( slice_realbox.hi(idim) > real_box.hi(idim) ) { + if ( slice_realbox.hi(idim) > real_box.hi(idim) ) { slice_realbox.setHi( idim, real_box.hi(idim)); - amrex::Print() << " slice hi is out of bounds." << - " Modified it in dimension " << idim << + amrex::Print() << " slice hi is out of bounds." << + " Modified it in dimension " << idim << " to be aligned with the domain box\n"; } - + // Factor to ensure index values computation depending on index type // double fac = ( 1.0 - SliceType[idim] )*dom_geom[0].CellSize(idim)*0.5; // if dimension is reduced to one cell length // @@ -264,58 +264,58 @@ CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, { slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) ); slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) ); - + if ( slice_cr_ratio[idim] > 1) slice_cr_ratio[idim] = 1; - + // check for interpolation -- compute index lo with floor and ceil - if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) >= fac ) { - slice_lo[idim] = floor( ( (slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) + fac ) ) + if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) >= fac ) { + slice_lo[idim] = floor( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac ) ) / dom_geom[0].CellSize(idim)) + fac * 1E-10); - slice_lo2[idim] = ceil( ( (slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) + fac) ) - / dom_geom[0].CellSize(idim)) - fac * 1E-10 ); - } - else { - slice_lo[idim] = round( (slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) ) ) + slice_lo2[idim] = ceil( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac) ) + / dom_geom[0].CellSize(idim)) - fac * 1E-10 ); + } + else { + slice_lo[idim] = round( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) / dom_geom[0].CellSize(idim)); - slice_lo2[idim] = ceil((slice_cc_nd_box.lo(idim) - - (real_box.lo(idim) ) ) + slice_lo2[idim] = ceil((slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) / dom_geom[0].CellSize(idim) ); } - + // flag for interpolation -- if reduced dimension location // // does not align with data point // - if ( slice_lo[idim] == slice_lo2[idim]) { + if ( slice_lo[idim] == slice_lo2[idim]) { if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) < fac ) { interp_lo[idim] = 1; } } - else { + else { interp_lo[idim] = 1; } - + // ncells = 1 if dimension is reduced // slice_hi[idim] = slice_lo[idim] + 1; } else { // moving realbox.lo and reabox.hi to nearest coarsenable grid point // - int index_lo = floor(((slice_realbox.lo(idim) + 1E-10 + int index_lo = floor(((slice_realbox.lo(idim) + 1E-10 - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); int index_hi = ceil(((slice_realbox.hi(idim) - 1E-10 - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); bool modify_cr = true; - - while ( modify_cr == true) { + + while ( modify_cr == true) { int lo_new = index_lo; - int hi_new = index_hi; + int hi_new = index_hi; int mod_lo = index_lo % slice_cr_ratio[idim]; int mod_hi = index_hi % slice_cr_ratio[idim]; modify_cr = false; - + // To ensure that the index.lo is coarsenable // if ( mod_lo > 0) { lo_new = index_lo - mod_lo; @@ -324,36 +324,36 @@ CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, if ( mod_hi > 0) { hi_new = index_hi + (slice_cr_ratio[idim] - mod_hi); } - - // If modified index.hi is > baselinebox.hi, move the point // + + // If modified index.hi is > baselinebox.hi, move the point // // to the previous coarsenable point // - if ( (hi_new * dom_geom[0].CellSize(idim)) + if ( (hi_new * dom_geom[0].CellSize(idim)) > real_box.hi(idim) - real_box.lo(idim) + dom_geom[0].CellSize(idim)*0.01 ) { hi_new = index_hi - mod_hi; } - - if ( (hi_new - lo_new) == 0 ){ + + if ( (hi_new - lo_new) == 0 ){ amrex::Print() << " Diagnostic Warning :: "; amrex::Print() << " Coarsening ratio "; amrex::Print() << slice_cr_ratio[idim] << " in dim "<< idim; amrex::Print() << "is leading to zero cells for slice."; amrex::Print() << " Thus reducing cr_ratio by half.\n"; - + slice_cr_ratio[idim] = slice_cr_ratio[idim]/2; modify_cr = true; } - - if ( modify_cr == false ) { + + if ( modify_cr == false ) { index_lo = lo_new; index_hi = hi_new; } slice_lo[idim] = index_lo; - slice_hi[idim] = index_hi - 1; // since default is cell-centered + slice_hi[idim] = index_hi - 1; // since default is cell-centered } - slice_realbox.setLo( idim, index_lo * dom_geom[0].CellSize(idim) + slice_realbox.setLo( idim, index_lo * dom_geom[0].CellSize(idim) + real_box.lo(idim) ); - slice_realbox.setHi( idim, index_hi * dom_geom[0].CellSize(idim) + slice_realbox.setHi( idim, index_hi * dom_geom[0].CellSize(idim) + real_box.lo(idim) ); slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) + fac ); slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) - fac ); @@ -363,14 +363,14 @@ CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, /* \brief - * This function is called if the coordinates of the slice do not align with data points - * \param interp_lo is an IntVect which is flagged as 1, if interpolation - is required in the dimension. + * This function is called if the coordinates of the slice do not align with data points + * \param interp_lo is an IntVect which is flagged as 1, if interpolation + is required in the dimension. */ -void +void InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, - Vector<Geometry> geom, int ncomp, int nghost, - IntVect slice_lo, IntVect slice_hi, IntVect SliceType, + Vector<Geometry> geom, int ncomp, int nghost, + IntVect slice_lo, IntVect slice_hi, IntVect SliceType, const RealBox real_box) { for (MFIter mfi(smf); mfi.isValid(); ++mfi) @@ -381,9 +381,9 @@ InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, const auto hi = amrex::ubound(bx); FArrayBox& fabox = smf[mfi]; - for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if ( interp_lo[idim] == 1 ) { - InterpolateLo( bx, fabox, slice_lo, geom, idim, SliceType, + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( interp_lo[idim] == 1 ) { + InterpolateLo( bx, fabox, slice_lo, geom, idim, SliceType, slice_realbox, 0, ncomp, nghost, real_box); } } @@ -391,10 +391,10 @@ InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, } -void -InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, - Vector<Geometry> geom, int idir, IntVect IndType, - RealBox slice_realbox, int srccomp, int ncomp, +void +InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, + Vector<Geometry> geom, int idir, IntVect IndType, + RealBox slice_realbox, int srccomp, int ncomp, int nghost, const RealBox real_box ) { auto fabarr = fabox.array(); @@ -428,7 +428,7 @@ InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, break; } case 1: - { + { if ( imin >= lo.y && imin <= lo.y) { for (int n = srccomp; n < srccomp+ncomp; ++n) { for (int k = lo.z; k <= hi.z; ++k) { diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index f79079b0e..6aa1f7a5b 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -750,7 +750,7 @@ WarpX::WriteJobInfo (const std::string& dir) const /* \brief * The raw slice data is written out in the plotfile format and can be visualized using yt. - * The slice data is written to diags/slice_plotfiles/pltXXXXX at the plotting interval. + * The slice data is written to diags/slice_plotfiles/pltXXXXX at the plotting interval. */ void WarpX::WriteSlicePlotFile () const @@ -763,36 +763,36 @@ WarpX::WriteSlicePlotFile () const VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(slice_plotfile_headerversion); rfs.emplace_back("raw_fields"); - + const int nlevels = finestLevel() + 1; - + // creating a temporary cell-centered dummy multifab // // to get around the issue of yt complaining about no field data // Vector< std::unique_ptr<MultiFab> > dummy_mf(nlevels); const DistributionMapping &dm2 = Efield_slice[0][0]->DistributionMap(); Vector<std::string> varnames; IntVect cc(AMREX_D_DECL(0,0,0)); - for (int lev = 0; lev < nlevels; ++lev) + for (int lev = 0; lev < nlevels; ++lev) { - dummy_mf[lev].reset(new MultiFab( - amrex::convert(Efield_slice[lev][0]->boxArray(),cc), + dummy_mf[lev].reset(new MultiFab( + amrex::convert(Efield_slice[lev][0]->boxArray(),cc), dm2, 1, 0 )); dummy_mf[lev]->setVal(0.0); - } + } amrex::WriteMultiLevelPlotfile(slice_plotfilename, nlevels, - GetVecOfConstPtrs(dummy_mf), + GetVecOfConstPtrs(dummy_mf), varnames, Geom(), t_new[0], istep, refRatio(), - "HyperCLaw-V1.1", + "HyperCLaw-V1.1", "Level_", "Cell", rfs); - for (int lev = 0; lev < nlevels; ++lev) + for (int lev = 0; lev < nlevels; ++lev) { const std::unique_ptr<MultiFab> empty_ptr; const std::string raw_spltname = slice_plotfilename + "/raw_fields"; amrex::Print() << " raw spltname " << raw_spltname << "\n"; const DistributionMapping &dm = Efield_slice[lev][0]->DistributionMap(); - + WriteRawField( *Efield_slice[lev][0], dm, raw_spltname, level_prefix, "Ex_slice", lev, 0); WriteRawField( *Efield_slice[lev][1], dm, raw_spltname, level_prefix, "Ey_slice", lev, 0); WriteRawField( *Efield_slice[lev][2], dm, raw_spltname, level_prefix, "Ez_slice", lev, 0); @@ -807,8 +807,8 @@ WarpX::WriteSlicePlotFile () const MultiFab rho_new(*rho_slice[lev], amrex::make_alias, 1, 1); WriteRawField( rho_new, dm, raw_spltname, level_prefix, "rho_slice", lev, 0); } - } - + } + WriteJobInfo(slice_plotfilename); WriteWarpXHeader(slice_plotfilename); @@ -817,7 +817,7 @@ WarpX::WriteSlicePlotFile () const } -void +void WarpX::InitializeSliceMultiFabs () { @@ -828,12 +828,12 @@ WarpX::InitializeSliceMultiFabs () current_slice.resize(nlevels); Efield_slice.resize(nlevels); Bfield_slice.resize(nlevels); - + } // To generate slice that inherits index type of underlying data // -void +void WarpX::SliceGenerationForDiagnostics () { @@ -841,20 +841,20 @@ WarpX::SliceGenerationForDiagnostics () dom_geom = Geom(); if (F_fp[0] ) { - F_slice[0] = CreateSlice( *F_fp[0].get(), dom_geom, slice_realbox, + F_slice[0] = CreateSlice( *F_fp[0].get(), dom_geom, slice_realbox, slice_cr_ratio ); } if (rho_fp[0]) { - rho_slice[0] = CreateSlice( *rho_fp[0].get(), dom_geom, slice_realbox, + rho_slice[0] = CreateSlice( *rho_fp[0].get(), dom_geom, slice_realbox, slice_cr_ratio ); } for (int idim = 0; idim < 3; ++idim) { - Efield_slice[0][idim] = CreateSlice( *Efield_fp[0][idim].get(), + Efield_slice[0][idim] = CreateSlice( *Efield_fp[0][idim].get(), dom_geom, slice_realbox, slice_cr_ratio ); - Bfield_slice[0][idim] = CreateSlice( *Bfield_fp[0][idim].get(), + Bfield_slice[0][idim] = CreateSlice( *Bfield_fp[0][idim].get(), dom_geom, slice_realbox, slice_cr_ratio ); - current_slice[0][idim] = CreateSlice( *current_fp[0][idim].get(), + current_slice[0][idim] = CreateSlice( *current_fp[0][idim].get(), dom_geom, slice_realbox, slice_cr_ratio ); } @@ -862,7 +862,7 @@ WarpX::SliceGenerationForDiagnostics () } -void +void WarpX::ClearSliceMultiFabs () { diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 8f800665d..a4ed803ee 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -148,7 +148,7 @@ WarpX::EvolveEM (int numsteps) } myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - + bool move_j = is_synchronized || to_make_plot || do_insitu; // If is_synchronized we need to shift j too so that next step we can evolve E by dt/2. // We might need to move j because we are going to make a plotfile. @@ -277,7 +277,7 @@ WarpX::EvolveEM (int numsteps) void WarpX::OneStep_nosub (Real cur_time) { - // Loop over species. For each ionizable species, create particles in + // Loop over species. For each ionizable species, create particles in // product species. mypc->doFieldIonization(); // Push particle from x^{n} to x^{n+1} @@ -349,7 +349,7 @@ WarpX::OneStep_sub1 (Real curtime) { // TODO: we could save some charge depositions - // Loop over species. For each ionizable species, create particles in + // Loop over species. For each ionizable species, create particles in // product species. mypc->doFieldIonization(); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H index a2511b6b7..50eb5c9b1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H @@ -17,7 +17,7 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm void InitializeSpectralCoefficients( const SpectralKSpace& spectral_kspace, - const amrex::DistributionMapping& dm, + const amrex::DistributionMapping& dm, const amrex::Real dt); // Redefine functions from base class diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 825d04dc2..bc7f90fac 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -20,7 +20,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm virtual int getRequiredNumberOfFields() const override final { return SpectralFieldIndex::n_fields; } - + void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const amrex::Real dt); diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 6fe5e3939..aee131324 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -149,7 +149,7 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, #else shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] ); #endif - + } } return shift_factor; diff --git a/Source/FieldSolver/openbc_poisson_solver.F90 b/Source/FieldSolver/openbc_poisson_solver.F90 index c18b1db24..8e962a2ef 100644 --- a/Source/FieldSolver/openbc_poisson_solver.F90 +++ b/Source/FieldSolver/openbc_poisson_solver.F90 @@ -4,7 +4,7 @@ module warpx_openbc_module implicit none integer, parameter :: idecomp = 0 ! 0=xyz, 1=xy, 2=yz, 3=xz, 4=x, 5=y, 6=z. - integer, parameter :: igfflag = 1 ! =0 for ordinary 1/r Green function; + integer, parameter :: igfflag = 1 ! =0 for ordinary 1/r Green function; ! =1 for integrated Green function integer, save :: gb_lo(3), gb_hi(3) diff --git a/Source/FieldSolver/solve_E_nodal.F90 b/Source/FieldSolver/solve_E_nodal.F90 index e5ac50b2a..4ed9bb845 100644 --- a/Source/FieldSolver/solve_E_nodal.F90 +++ b/Source/FieldSolver/solve_E_nodal.F90 @@ -8,7 +8,7 @@ module warpx_ES_solve_E_nodal contains ! This routine computes the node-centered electric field given a node-centered phi. -! The gradient is computed using 2nd-order centered differences. It assumes the +! The gradient is computed using 2nd-order centered differences. It assumes the ! Boundary conditions have already been set and that you have two rows of ghost cells ! for phi and one row of ghost cells for Ex, Ey, and Ez. ! Note that this routine includes the minus sign in E = - grad phi. @@ -35,7 +35,7 @@ contains do k = lo(3)-1, hi(3)+1 do j = lo(2)-1, hi(2)+1 do i = lo(1)-1, hi(1)+1 - + Ex(i,j,k) = fac(1) * (phi(i-1,j,k) - phi(i+1,j,k)) Ey(i,j,k) = fac(2) * (phi(i,j-1,k) - phi(i,j+1,k)) Ez(i,j,k) = fac(3) * (phi(i,j,k-1) - phi(i,j,k+1)) @@ -61,7 +61,7 @@ contains do j = lo(2)-1, hi(2)+1 do i = lo(1)-1, hi(1)+1 - + Ex(i,j) = fac(1) * (phi(i-1,j) - phi(i+1,j)) Ey(i,j) = fac(2) * (phi(i,j-1) - phi(i,j+1)) diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index 4fcfec14b..150ca8e08 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -8,7 +8,7 @@ class BilinearFilter : public Filter public: BilinearFilter() = default; - + void ComputeStencils(); amrex::IntVect npass_each_dir; diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index 68ebde687..23abaf8bc 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -23,8 +23,8 @@ namespace { new_s[0] = 0.5 * old_s[0]; if (1<jmax) new_s[0] += 0.5 * old_s[1]; loc = 0.; - // For each element j, apply the filter to - // old_s to get new_s[j]. loc stores the tmp + // For each element j, apply the filter to + // old_s to get new_s[j]. loc stores the tmp // filtered value. for(int j=1; j<jmax+1; j++){ loc = 0.5 * old_s[j]; diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index eaa0498c9..5eb84b96a 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -11,13 +11,13 @@ public: // Apply stencil on MultiFab. // Guard cells are handled inside this function - void ApplyStencil(amrex::MultiFab& dstmf, - const amrex::MultiFab& srcmf, int scomp=0, + void ApplyStencil(amrex::MultiFab& dstmf, + const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); // Apply stencil on a FabArray. - void ApplyStencil (amrex::FArrayBox& dstfab, - const amrex::FArrayBox& srcfab, const amrex::Box& tbx, + void ApplyStencil (amrex::FArrayBox& dstfab, + const amrex::FArrayBox& srcfab, const amrex::Box& tbx, int scomp=0, int dcomp=0, int ncomp=10000); // public for cuda diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index f8a2dd6c2..5d3c14c14 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -59,7 +59,7 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); @@ -106,11 +106,11 @@ void Filter::DoFilter (const Box& tbx, for (int iz=0; iz < slen_local.z; ++iz){ for (int iy=0; iy < slen_local.y; ++iy){ for (int ix=0; ix < slen_local.x; ++ix){ -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) Real sss = sx[ix]*sy[iy]*sz[iz]; #else Real sss = sx[ix]*sz[iy]; -#endif +#endif #if (AMREX_SPACEDIM == 3) d += sss*( tmp(i-ix,j-iy,k-iz,scomp+n) +tmp(i+ix,j-iy,k-iz,scomp+n) @@ -179,7 +179,7 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); @@ -220,7 +220,7 @@ void Filter::DoFilter (const Box& tbx, for (int iz=0; iz < slen.z; ++iz){ for (int iy=0; iy < slen.y; ++iy){ for (int ix=0; ix < slen.x; ++ix){ -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) Real sss = sx[ix]*sy[iy]*sz[iz]; #else Real sss = sx[ix]*sz[iy]; diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H index a53039dfa..4116c3a48 100644 --- a/Source/Filter/NCIGodfreyFilter.H +++ b/Source/Filter/NCIGodfreyFilter.H @@ -20,7 +20,7 @@ public: static constexpr int stencil_width = 4; private: - + godfrey_coeff_set coeff_set; amrex::Real cdtodz; int l_lower_order_in_v; diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 3f589260a..25b5e67fd 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -20,7 +20,7 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdt #else stencil_length_each_dir = {1,5}; slen = {1,5,1}; -#endif +#endif } void NCIGodfreyFilter::ComputeStencils(){ @@ -28,12 +28,12 @@ void NCIGodfreyFilter::ComputeStencils(){ #if ( AMREX_SPACEDIM == 3 ) slen.z==5, #else - slen.y==5, + slen.y==5, #endif "ERROR: NCI filter requires 5 points in z"); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_lower_order_in_v==1, + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_lower_order_in_v==1, "ERROR: NCI corrector requires l_lower_order_in_v=1, i.e., Galerkin scheme"); - + // Interpolate coefficients from the table, and store into prestencil. int index = tab_length*cdtodz; index = min(index, tab_length-2); @@ -42,10 +42,10 @@ void NCIGodfreyFilter::ComputeStencils(){ Real prestencil[4]; for(int i=0; i<tab_width; i++){ if (coeff_set == godfrey_coeff_set::Ex_Ey_Bz) { - prestencil[i] = (1-weight_right)*table_nci_godfrey_Ex_Ey_Bz[index ][i] + + prestencil[i] = (1-weight_right)*table_nci_godfrey_Ex_Ey_Bz[index ][i] + weight_right*table_nci_godfrey_Ex_Ey_Bz[index+1][i]; } else if (coeff_set == godfrey_coeff_set::Bx_By_Ez) { - prestencil[i] = (1-weight_right)*table_nci_godfrey_Bx_By_Ez[index ][i] + + prestencil[i] = (1-weight_right)*table_nci_godfrey_Bx_By_Ez[index ][i] + weight_right*table_nci_godfrey_Bx_By_Ez[index+1][i]; } } @@ -58,7 +58,7 @@ void NCIGodfreyFilter::ComputeStencils(){ stencil_z[3] = -( 4*prestencil[2] + 8*prestencil[3]) / 256; stencil_z[4] = ( 1*prestencil[3]) / 256; - // Compute stencil_x and stencil_y (no filter in these directions, + // Compute stencil_x and stencil_y (no filter in these directions, // so only 1 coeff, equal to 1) stencil_x.resize(1); stencil_x[0] = 1.; @@ -67,12 +67,12 @@ void NCIGodfreyFilter::ComputeStencils(){ stencil_y[0] = 1.; #endif - // Due to the way Filter::DoFilter() is written, + // Due to the way Filter::DoFilter() is written, // coefficient 0 has to be /2 stencil_x[0] /= 2.; #if (AMREX_SPACEDIM == 3) stencil_y[0] /= 2.; #endif stencil_z[0] /= 2.; - + } diff --git a/Source/Initialization/CustomDensityProb.H b/Source/Initialization/CustomDensityProb.H index b00830e6c..c5c159ee8 100644 --- a/Source/Initialization/CustomDensityProb.H +++ b/Source/Initialization/CustomDensityProb.H @@ -27,7 +27,7 @@ struct InjectorDensityCustom } } - // Return density at given position, using user-defined parameters + // Return density at given position, using user-defined parameters // stored in p. AMREX_GPU_HOST_DEVICE amrex::Real diff --git a/Source/Initialization/InjectorDensity.H b/Source/Initialization/InjectorDensity.H index afcc23c4b..de5f9cb4c 100644 --- a/Source/Initialization/InjectorDensity.H +++ b/Source/Initialization/InjectorDensity.H @@ -95,9 +95,9 @@ private: amrex::Real* p; }; -// Base struct for density injector. -// InjectorDensity contains a union (called Object) that holds any one -// instance of: +// Base struct for density injector. +// InjectorDensity contains a union (called Object) that holds any one +// instance of: // - InjectorDensityConstant : to generate constant density; // - InjectorDensityParser : to generate density from parser; // - InjectorDensityCustom : to generate density from custom profile; diff --git a/Source/Initialization/InjectorDensity.cpp b/Source/Initialization/InjectorDensity.cpp index 21af95f04..ddf6d1a18 100644 --- a/Source/Initialization/InjectorDensity.cpp +++ b/Source/Initialization/InjectorDensity.cpp @@ -49,7 +49,7 @@ InjectorDensityPredefined::InjectorDensityPredefined ( ParmParse pp(a_species_name); std::vector<amrex::Real> v; - // Read parameters for the predefined plasma profile, + // Read parameters for the predefined plasma profile, // and store them in managed memory pp.getarr("predefined_profile_params", v); p = static_cast<amrex::Real*> diff --git a/Source/Initialization/InjectorMomentum.H b/Source/Initialization/InjectorMomentum.H index 7dced1336..ba0d26fc4 100644 --- a/Source/Initialization/InjectorMomentum.H +++ b/Source/Initialization/InjectorMomentum.H @@ -23,7 +23,7 @@ private: amrex::Real m_ux, m_uy, m_uz; }; -// struct whose getMomentum returns momentum for 1 particle, from random +// struct whose getMomentum returns momentum for 1 particle, from random // gaussian distribution. struct InjectorMomentumGaussian { @@ -85,9 +85,9 @@ struct InjectorMomentumParser GpuParser m_ux_parser, m_uy_parser, m_uz_parser; }; -// Base struct for momentum injector. -// InjectorMomentum contains a union (called Object) that holds any one -// instance of: +// Base struct for momentum injector. +// InjectorMomentum contains a union (called Object) that holds any one +// instance of: // - InjectorMomentumConstant : to generate constant density; // - InjectorMomentumGaussian : to generate gaussian distribution; // - InjectorMomentumRadialExpansion: to generate radial expansion; @@ -190,7 +190,7 @@ private: Type type; // An instance of union Object constructs and stores any one of - // the objects declared (constant or custom or gaussian or + // the objects declared (constant or custom or gaussian or // radial_expansion or parser). union Object { Object (InjectorMomentumConstant*, diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index 0e43f889e..f8f16746c 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -25,7 +25,7 @@ struct InjectorPositionRegular // i_part: particle number within the cell, required to evenly space // particles within the cell. - // ref_fac: the number of particles evenly-spaced within a cell + // ref_fac: the number of particles evenly-spaced within a cell // is a_ppc*(ref_fac**AMREX_SPACEDIM). AMREX_GPU_HOST_DEVICE amrex::XDim3 @@ -47,9 +47,9 @@ private: amrex::Dim3 ppc; }; -// Base struct for position injector. -// InjectorPosition contains a union (called Object) that holds any one -// instance of: +// Base struct for position injector. +// InjectorPosition contains a union (called Object) that holds any one +// instance of: // - InjectorPositionRandom : to generate random distribution; // - InjectorPositionRegular: to generate regular distribution. // The choice is made at runtime, depending in the constructor called. diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index 875fdd4ad..a944165d6 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -96,9 +96,9 @@ protected: std::unique_ptr<InjectorPosition> inj_pos; std::unique_ptr<InjectorDensity > inj_rho; std::unique_ptr<InjectorMomentum> inj_mom; - + void parseDensity (amrex::ParmParse& pp); - void parseMomentum (amrex::ParmParse& pp); + void parseMomentum (amrex::ParmParse& pp); }; #endif diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 2fd42f168..782282a26 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -18,12 +18,12 @@ public: virtual void InitData () final; -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, amrex::Real t, amrex::Real dt) { BL_ASSERT(false); } #endif // WARPX_DO_ELECTROSTATIC - + virtual void Evolve (int lev, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, const amrex::MultiFab&, @@ -64,7 +64,7 @@ public: const amrex::Real dt, const int thread_num); protected: - + std::string laser_name; private: @@ -107,11 +107,11 @@ private: // Theoretical position of the antenna. Used if do_continuous_injection=1. // Track the position of the antenna until it enters the simulation domain. amrex::Vector<amrex::Real> updated_position; - + void ComputeSpacing (int lev, amrex::Real& Sx, amrex::Real& Sy) const; void ComputeWeightMobility (amrex::Real Sx, amrex::Real Sy); void InitData (int lev); - // Inject the laser antenna during the simulation, if it started + // Inject the laser antenna during the simulation, if it started // outside of the simulation domain and enters it. void ContinuousInjection(const amrex::RealBox& injection_box) override; // Update position of the antenna diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 9699ee06b..d0ff51313 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -27,7 +27,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, charge = 1.0; mass = std::numeric_limits<Real>::max(); do_boosted_frame_diags = 0; - + ParmParse pp(laser_name); // Parse the type of laser profile and set the corresponding flag `profile` @@ -166,7 +166,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, // If laser antenna initially outside of the box, store its theoretical // position in z_antenna_th updated_position = position; - + // Sanity checks int dir = WarpX::moving_window_dir; std::vector<Real> windir(3, 0.0); @@ -175,16 +175,16 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, #else windir[dir] = 1.0; #endif - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - (nvec[0]-windir[0]) + (nvec[1]-windir[1]) + (nvec[2]-windir[2]) + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (nvec[0]-windir[0]) + (nvec[1]-windir[1]) + (nvec[2]-windir[2]) < 1.e-12, "do_continous_injection for laser particle only works" + " if moving window direction and laser propagation direction are the same"); if ( WarpX::gamma_boost>1 ){ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + - (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + + (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + + (WarpX::boost_direction[1]-0)*(WarpX::boost_direction[1]-0) + (WarpX::boost_direction[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, - "do_continous_injection for laser particle only works if " + + "do_continous_injection for laser particle only works if " + "warpx.boost_direction = z. TODO: all directions."); } } @@ -192,13 +192,13 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, /* \brief Check if laser particles enter the box, and inject if necessary. * \param injection_box: a RealBox where particles should be injected. - */ + */ void LaserParticleContainer::ContinuousInjection (const RealBox& injection_box) { // Input parameter injection_box contains small box where injection // should occur. - // So far, LaserParticleContainer::laser_injection_box contains the + // So far, LaserParticleContainer::laser_injection_box contains the // outdated full problem domain at t=0. // Convert updated_position to Real* to use RealBox::contains(). @@ -258,7 +258,7 @@ LaserParticleContainer::InitData (int lev) ComputeSpacing(lev, S_X, S_Y); ComputeWeightMobility(S_X, S_Y); - // LaserParticleContainer::position contains the initial position of the + // LaserParticleContainer::position contains the initial position of the // laser antenna. In the boosted frame, the antenna is moving. // Update its position with updated_position. if (do_continuous_injection){ @@ -459,7 +459,7 @@ LaserParticleContainer::Evolve (int lev, if (rho) { int* AMREX_RESTRICT ion_lev = nullptr; - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (crho) { DepositCharge(pti, wp, ion_lev, crho, 0, np_current, @@ -606,7 +606,7 @@ LaserParticleContainer::PushP (int lev, Real dt, * \param np: number of laser particles * \param thread_num: thread number * \param pplane_Xp, pplane_Yp: pointers to arrays of particle positions - * in laser plane coordinate. + * in laser plane coordinate. */ void LaserParticleContainer::calculate_laser_plane_coordinates ( @@ -633,16 +633,16 @@ LaserParticleContainer::calculate_laser_plane_coordinates ( np, [=] AMREX_GPU_DEVICE (int i) { #if (AMREX_SPACEDIM == 3) - pplane_Xp[i] = + pplane_Xp[i] = tmp_u_X_0 * (xp[i] - tmp_position_0) + tmp_u_X_1 * (yp[i] - tmp_position_1) + tmp_u_X_2 * (zp[i] - tmp_position_2); - pplane_Yp[i] = + pplane_Yp[i] = tmp_u_Y_0 * (xp[i] - tmp_position_0) + tmp_u_Y_1 * (yp[i] - tmp_position_1) + tmp_u_Y_2 * (zp[i] - tmp_position_2); #elif (AMREX_SPACEDIM == 2) - pplane_Xp[i] = + pplane_Xp[i] = tmp_u_X_0 * (xp[i] - tmp_position_0) + tmp_u_X_2 * (zp[i] - tmp_position_2); pplane_Yp[i] = 0.; @@ -677,7 +677,7 @@ LaserParticleContainer::update_laser_particle( Real tmp_nvec_0 = nvec[0]; Real tmp_nvec_1 = nvec[1]; Real tmp_nvec_2 = nvec[2]; - + // Copy member variables to tmp copies for GPU runs. Real tmp_mobility = mobility; Real gamma_boost = WarpX::gamma_boost; diff --git a/Source/Laser/LaserProfiles.cpp b/Source/Laser/LaserProfiles.cpp index 704107c2d..69804b17c 100644 --- a/Source/Laser/LaserProfiles.cpp +++ b/Source/Laser/LaserProfiles.cpp @@ -9,14 +9,14 @@ using namespace amrex; * Both Xp and Yp are given in laser plane coordinate. * For each particle with position Xp and Yp, this routine computes the * amplitude of the laser electric field, stored in array amplitude. - * + * * \param np: number of laser particles * \param Xp: pointer to first component of positions of laser particles * \param Yp: pointer to second component of positions of laser particles * \param t: Current physical time * \param amplitude: pointer to array of field amplitude. */ -void +void LaserParticleContainer::gaussian_laser_profile ( const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, Real t, Real * AMREX_RESTRICT const amplitude) @@ -34,9 +34,9 @@ LaserParticleContainer::gaussian_laser_profile ( // Time stretching due to STCs and phi2 complex envelope // (1 if zeta=0, beta=0, phi2=0) - const Complex stretch_factor = 1. + 4. * + const Complex stretch_factor = 1. + 4. * (zeta+beta*profile_focal_distance) * (zeta+beta*profile_focal_distance) - * (inv_tau2*inv_complex_waist_2) + + * (inv_tau2*inv_complex_waist_2) + 2.*I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; // Amplitude and monochromatic oscillations @@ -59,10 +59,10 @@ LaserParticleContainer::gaussian_laser_profile ( Real tmp_profile_focal_distance = profile_focal_distance; // Loop through the macroparticle to calculate the proper amplitude amrex::ParallelFor( - np, + np, [=] AMREX_GPU_DEVICE (int i) { const Complex stc_exponent = 1./stretch_factor * inv_tau2 * - MathFunc::pow((t - tmp_profile_t_peak - + MathFunc::pow((t - tmp_profile_t_peak - tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - 2.*I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) *( tmp_zeta - tmp_beta*tmp_profile_focal_distance ) * inv_complex_waist_2),2); @@ -81,14 +81,14 @@ LaserParticleContainer::gaussian_laser_profile ( * Both Xp and Yp are given in laser plane coordinate. * For each particle with position Xp and Yp, this routine computes the * amplitude of the laser electric field, stored in array amplitude. - * + * * \param np: number of laser particles * \param Xp: pointer to first component of positions of laser particles * \param Yp: pointer to second component of positions of laser particles * \param t: Current physical time * \param amplitude: pointer to array of field amplitude. */ -void +void LaserParticleContainer::harris_laser_profile ( const int np, Real const * AMREX_RESTRICT const Xp, Real const * AMREX_RESTRICT const Yp, Real t, Real * AMREX_RESTRICT const amplitude) @@ -96,14 +96,14 @@ LaserParticleContainer::harris_laser_profile ( // This function uses the Harris function as the temporal profile of the pulse const Real omega0 = 2.*MathConst::pi*PhysConst::c/wavelength; const Real zR = MathConst::pi * profile_waist*profile_waist / wavelength; - const Real wz = profile_waist * + const Real wz = profile_waist * std::sqrt(1. + profile_focal_distance*profile_focal_distance/zR*zR); const Real inv_wz_2 = 1./(wz*wz); Real inv_Rz; - if (profile_focal_distance == 0.){ + if (profile_focal_distance == 0.){ inv_Rz = 0.; } else { - inv_Rz = -profile_focal_distance / + inv_Rz = -profile_focal_distance / ( profile_focal_distance*profile_focal_distance + zR*zR ); } const Real arg_env = 2.*MathConst::pi*t/profile_duration; @@ -111,8 +111,8 @@ LaserParticleContainer::harris_laser_profile ( // time envelope is given by the Harris function Real time_envelope = 0.; if (t < profile_duration) - time_envelope = 1./32. * (10. - 15.*std::cos(arg_env) + - 6.*std::cos(2.*arg_env) - + time_envelope = 1./32. * (10. - 15.*std::cos(arg_env) + + 6.*std::cos(2.*arg_env) - std::cos(3.*arg_env)); // Copy member variables to tmp copies for GPU runs. @@ -121,7 +121,7 @@ LaserParticleContainer::harris_laser_profile ( amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) { - const Real space_envelope = + const Real space_envelope = std::exp(- ( Xp[i]*Xp[i] + Yp[i]*Yp[i] ) * inv_wz_2); const Real arg_osc = omega0*t - omega0/PhysConst::c* (Xp[i]*Xp[i] + Yp[i]*Yp[i]) * inv_Rz / 2.; diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 9b3baafe1..c737591c7 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -33,11 +33,11 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa { if (ParallelDescriptor::NProcs() == 1) return; -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC AMREX_ALWAYS_ASSERT(masks[lev] == nullptr); AMREX_ALWAYS_ASSERT(gather_masks[lev] == nullptr); #endif // WARPX_DO_ELECTROSTATIC - + // Fine patch const auto& period = Geom(lev).periodicity(); @@ -203,7 +203,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa auto pmf = std::unique_ptr<iMultiFab>(new iMultiFab(current_buffer_masks[lev]->boxArray(), dm, current_buffer_masks[lev]->nComp(), ng)); // pmf->ParallelCopy(*current_buffer_masks[lev], 0, 0, current_buffer_masks[lev]->nComp(), ng, ng); - current_buffer_masks[lev] = std::move(pmf); + current_buffer_masks[lev] = std::move(pmf); } if (gather_buffer_masks[lev]) { @@ -211,7 +211,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa auto pmf = std::unique_ptr<iMultiFab>(new iMultiFab(gather_buffer_masks[lev]->boxArray(), dm, gather_buffer_masks[lev]->nComp(), ng)); // pmf->ParallelCopy(*gather_buffer_masks[lev], 0, 0, gather_buffer_masks[lev]->nComp(), ng, ng); - gather_buffer_masks[lev] = std::move(pmf); + gather_buffer_masks[lev] = std::move(pmf); } } diff --git a/Source/Parser/GpuParser.H b/Source/Parser/GpuParser.H index 1533ee6b9..e49671e06 100644 --- a/Source/Parser/GpuParser.H +++ b/Source/Parser/GpuParser.H @@ -5,7 +5,7 @@ #include <AMReX_Gpu.H> // When compiled for CPU, wrap WarpXParser and enable threading. -// When compiled for GPU, store one copy of the parser in +// When compiled for GPU, store one copy of the parser in // CUDA managed memory for __device__ code, and one copy of the parser // in CUDA managed memory for __host__ code. This way, the parser can be // efficiently called from both host and device. diff --git a/Source/Parser/GpuParser.cpp b/Source/Parser/GpuParser.cpp index db1c2287d..5078b498b 100644 --- a/Source/Parser/GpuParser.cpp +++ b/Source/Parser/GpuParser.cpp @@ -28,7 +28,7 @@ GpuParser::GpuParser (WarpXParser const& wp) wp_parser_regvar(&m_cpu_parser, "x", &(m_var.x)); wp_parser_regvar(&m_cpu_parser, "y", &(m_var.y)); wp_parser_regvar(&m_cpu_parser, "z", &(m_var.z)); - + #else // not defined AMREX_USE_GPU #ifdef _OPENMP diff --git a/Source/Parser/wp_parser.lex.c b/Source/Parser/wp_parser.lex.c index 0a9c58ee9..09a1b2db8 100644 --- a/Source/Parser/wp_parser.lex.c +++ b/Source/Parser/wp_parser.lex.c @@ -34,7 +34,7 @@ #if defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L /* C99 says to define __STDC_LIMIT_MACROS before including stdint.h, - * if you want the limit (max/min) macros for int types. + * if you want the limit (max/min) macros for int types. */ #ifndef __STDC_LIMIT_MACROS #define __STDC_LIMIT_MACROS 1 @@ -51,7 +51,7 @@ typedef uint32_t flex_uint32_t; typedef signed char flex_int8_t; typedef short int flex_int16_t; typedef int flex_int32_t; -typedef unsigned char flex_uint8_t; +typedef unsigned char flex_uint8_t; typedef unsigned short int flex_uint16_t; typedef unsigned int flex_uint32_t; @@ -162,10 +162,10 @@ extern FILE *yyin, *yyout; #define EOB_ACT_CONTINUE_SCAN 0 #define EOB_ACT_END_OF_FILE 1 #define EOB_ACT_LAST_MATCH 2 - + #define YY_LESS_LINENO(n) #define YY_LINENO_REWIND_TO(ptr) - + /* Return all but the first "n" matched characters back to the input stream. */ #define yyless(n) \ do \ @@ -590,9 +590,9 @@ extern int yywrap ( void ); #endif #ifndef YY_NO_UNPUT - + static void yyunput ( int c, char *buf_ptr ); - + #endif #ifndef yytext_ptr @@ -719,7 +719,7 @@ YY_DECL yy_state_type yy_current_state; char *yy_cp, *yy_bp; int yy_act; - + if ( !(yy_init) ) { (yy_init) = 1; @@ -1300,7 +1300,7 @@ static int yy_get_next_buffer (void) { yy_state_type yy_current_state; char *yy_cp; - + yy_current_state = (yy_start); for ( yy_cp = (yytext_ptr) + YY_MORE_ADJ; yy_cp < (yy_c_buf_p); ++yy_cp ) @@ -1356,7 +1356,7 @@ static int yy_get_next_buffer (void) static void yyunput (int c, char * yy_bp ) { char *yy_cp; - + yy_cp = (yy_c_buf_p); /* undo effects of setting up yytext */ @@ -1401,7 +1401,7 @@ static int yy_get_next_buffer (void) { int c; - + *(yy_c_buf_p) = (yy_hold_char); if ( *(yy_c_buf_p) == YY_END_OF_BUFFER_CHAR ) @@ -1468,12 +1468,12 @@ static int yy_get_next_buffer (void) /** Immediately switch to a different input stream. * @param input_file A readable stream. - * + * * @note This function does not reset the start condition to @c INITIAL . */ void yyrestart (FILE * input_file ) { - + if ( ! YY_CURRENT_BUFFER ){ yyensure_buffer_stack (); YY_CURRENT_BUFFER_LVALUE = @@ -1486,11 +1486,11 @@ static int yy_get_next_buffer (void) /** Switch to a different input buffer. * @param new_buffer The new input buffer. - * + * */ void yy_switch_to_buffer (YY_BUFFER_STATE new_buffer ) { - + /* TODO. We should be able to replace this entire function body * with * yypop_buffer_state(); @@ -1530,13 +1530,13 @@ static void yy_load_buffer_state (void) /** Allocate and initialize an input buffer state. * @param file A readable stream. * @param size The character buffer size in bytes. When in doubt, use @c YY_BUF_SIZE. - * + * * @return the allocated buffer state. */ YY_BUFFER_STATE yy_create_buffer (FILE * file, int size ) { YY_BUFFER_STATE b; - + b = (YY_BUFFER_STATE) yyalloc( sizeof( struct yy_buffer_state ) ); if ( ! b ) YY_FATAL_ERROR( "out of dynamic memory in yy_create_buffer()" ); @@ -1559,11 +1559,11 @@ static void yy_load_buffer_state (void) /** Destroy the buffer. * @param b a buffer created with yy_create_buffer() - * + * */ void yy_delete_buffer (YY_BUFFER_STATE b ) { - + if ( ! b ) return; @@ -1584,7 +1584,7 @@ static void yy_load_buffer_state (void) { int oerrno = errno; - + yy_flush_buffer( b ); b->yy_input_file = file; @@ -1600,13 +1600,13 @@ static void yy_load_buffer_state (void) } b->yy_is_interactive = file ? (isatty( fileno(file) ) > 0) : 0; - + errno = oerrno; } /** Discard all buffered characters. On the next scan, YY_INPUT will be called. * @param b the buffer state to be flushed, usually @c YY_CURRENT_BUFFER. - * + * */ void yy_flush_buffer (YY_BUFFER_STATE b ) { @@ -1635,7 +1635,7 @@ static void yy_load_buffer_state (void) * the current state. This function will allocate the stack * if necessary. * @param new_buffer The new state. - * + * */ void yypush_buffer_state (YY_BUFFER_STATE new_buffer ) { @@ -1665,7 +1665,7 @@ void yypush_buffer_state (YY_BUFFER_STATE new_buffer ) /** Removes and deletes the top of the stack, if present. * The next element becomes the new top. - * + * */ void yypop_buffer_state (void) { @@ -1689,7 +1689,7 @@ void yypop_buffer_state (void) static void yyensure_buffer_stack (void) { yy_size_t num_to_alloc; - + if (!(yy_buffer_stack)) { /* First allocation is just for 2 elements, since we don't know if this @@ -1732,13 +1732,13 @@ static void yyensure_buffer_stack (void) /** Setup the input buffer state to scan directly from a user-specified character buffer. * @param base the character buffer * @param size the size in bytes of the character buffer - * + * * @return the newly allocated buffer state object. */ YY_BUFFER_STATE yy_scan_buffer (char * base, yy_size_t size ) { YY_BUFFER_STATE b; - + if ( size < 2 || base[size-2] != YY_END_OF_BUFFER_CHAR || base[size-1] != YY_END_OF_BUFFER_CHAR ) @@ -1767,14 +1767,14 @@ YY_BUFFER_STATE yy_scan_buffer (char * base, yy_size_t size ) /** Setup the input buffer state to scan a string. The next call to yylex() will * scan from a @e copy of @a str. * @param yystr a NUL-terminated string to scan - * + * * @return the newly allocated buffer state object. * @note If you want to scan bytes that may contain NUL values, then use * yy_scan_bytes() instead. */ YY_BUFFER_STATE yy_scan_string (const char * yystr ) { - + return yy_scan_bytes( yystr, (int) strlen(yystr) ); } @@ -1782,7 +1782,7 @@ YY_BUFFER_STATE yy_scan_string (const char * yystr ) * scan from a @e copy of @a bytes. * @param yybytes the byte buffer to scan * @param _yybytes_len the number of bytes in the buffer pointed to by @a bytes. - * + * * @return the newly allocated buffer state object. */ YY_BUFFER_STATE yy_scan_bytes (const char * yybytes, int _yybytes_len ) @@ -1791,7 +1791,7 @@ YY_BUFFER_STATE yy_scan_bytes (const char * yybytes, int _yybytes_len ) char *buf; yy_size_t n; int i; - + /* Get memory for full buffer, including space for trailing EOB's. */ n = (yy_size_t) (_yybytes_len + 2); buf = (char *) yyalloc( n ); @@ -1845,16 +1845,16 @@ static void yynoreturn yy_fatal_error (const char* msg ) /* Accessor methods (get/set functions) to struct members. */ /** Get the current line number. - * + * */ int yyget_lineno (void) { - + return yylineno; } /** Get the input stream. - * + * */ FILE *yyget_in (void) { @@ -1862,7 +1862,7 @@ FILE *yyget_in (void) } /** Get the output stream. - * + * */ FILE *yyget_out (void) { @@ -1870,7 +1870,7 @@ FILE *yyget_out (void) } /** Get the length of the current token. - * + * */ int yyget_leng (void) { @@ -1878,7 +1878,7 @@ int yyget_leng (void) } /** Get the current token. - * + * */ char *yyget_text (void) @@ -1888,18 +1888,18 @@ char *yyget_text (void) /** Set the current line number. * @param _line_number line number - * + * */ void yyset_lineno (int _line_number ) { - + yylineno = _line_number; } /** Set the input stream. This does not discard the current * input buffer. * @param _in_str A readable stream. - * + * * @see yy_switch_to_buffer */ void yyset_in (FILE * _in_str ) @@ -1953,7 +1953,7 @@ static int yy_init_globals (void) /* yylex_destroy is for both reentrant and non-reentrant scanners. */ int yylex_destroy (void) { - + /* Pop the buffer stack, destroying each element. */ while(YY_CURRENT_BUFFER){ yy_delete_buffer( YY_CURRENT_BUFFER ); @@ -1979,7 +1979,7 @@ int yylex_destroy (void) #ifndef yytext_ptr static void yy_flex_strncpy (char* s1, const char * s2, int n ) { - + int i; for ( i = 0; i < n; ++i ) s1[i] = s2[i]; @@ -2004,7 +2004,7 @@ void *yyalloc (yy_size_t size ) void *yyrealloc (void * ptr, yy_size_t size ) { - + /* The cast to (char *) in the following accommodates both * implementations that use char* generic pointers, and those * that use void* generic pointers. It works with the latter diff --git a/Source/Parser/wp_parser.lex.h b/Source/Parser/wp_parser.lex.h index 4fd613978..557add970 100644 --- a/Source/Parser/wp_parser.lex.h +++ b/Source/Parser/wp_parser.lex.h @@ -38,7 +38,7 @@ #if defined (__STDC_VERSION__) && __STDC_VERSION__ >= 199901L /* C99 says to define __STDC_LIMIT_MACROS before including stdint.h, - * if you want the limit (max/min) macros for int types. + * if you want the limit (max/min) macros for int types. */ #ifndef __STDC_LIMIT_MACROS #define __STDC_LIMIT_MACROS 1 @@ -55,7 +55,7 @@ typedef uint32_t flex_uint32_t; typedef signed char flex_int8_t; typedef short int flex_int16_t; typedef int flex_int32_t; -typedef unsigned char flex_uint8_t; +typedef unsigned char flex_uint8_t; typedef unsigned short int flex_uint16_t; typedef unsigned int flex_uint32_t; diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index f02eb1033..b9210e67c 100755 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -18,8 +18,8 @@ * /param q : species charge. */ template <int depos_order> -void doChargeDepositionShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, +void doChargeDepositionShapeN(const amrex::Real * const xp, + const amrex::Real * const yp, const amrex::Real * const zp, const amrex::Real * const wp, const int * const ion_lev, @@ -69,7 +69,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, amrex::Real sx[depos_order + 1]; // i: leftmost grid point (node-centered) that the particle touches const int i = compute_shape_factor<depos_order>(sx, x); - + #if (defined WARPX_DIM_3D) // y direction const amrex::Real y = (yp[ip] - ymin)*dyi; @@ -86,7 +86,7 @@ void doChargeDepositionShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( - &rho_arr(lo.x+i+ix, lo.y+k+iz, 0), + &rho_arr(lo.x+i+ix, lo.y+k+iz, 0), sx[ix]*sz[iz]*wq); } } diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index c7dfde75a..7a96dab9a 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -24,8 +24,8 @@ * /param q : species charge. */ template <int depos_order> -void doDepositionShapeN(const amrex::Real * const xp, - const amrex::Real * const yp, +void doDepositionShapeN(const amrex::Real * const xp, + const amrex::Real * const yp, const amrex::Real * const zp, const amrex::Real * const wp, const amrex::Real * const uxp, @@ -39,7 +39,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const std::array<amrex::Real,3>& dx, const std::array<amrex::Real, 3> xyzmin, const amrex::Dim3 lo, - const amrex::Real stagger_shift, + const amrex::Real stagger_shift, const amrex::Real q) { // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer @@ -77,7 +77,7 @@ void doDepositionShapeN(const amrex::Real * const xp, const amrex::Real vx = uxp[ip]*gaminv; const amrex::Real vy = uyp[ip]*gaminv; const amrex::Real vz = uzp[ip]*gaminv; - // wqx, wqy wqz are particle current in each direction + // wqx, wqy wqz are particle current in each direction #if (defined WARPX_DIM_RZ) // In RZ, wqx is actually wqr, and wqy is wqtheta // Convert to cylinderical at the mid point @@ -117,7 +117,7 @@ void doDepositionShapeN(const amrex::Real * const xp, amrex::Real sx0[depos_order + 1]; // j0: leftmost grid point (cell-centered) that the particle touches const int j0 = compute_shape_factor<depos_order>(sx0, xmid-stagger_shift); - + #if (defined WARPX_DIM_3D) // y direction const amrex::Real ymid= (yp[ip]-ymin)*dyi-dts2dy*vy; @@ -138,13 +138,13 @@ void doDepositionShapeN(const amrex::Real * const xp, for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::Add( - &jx_arr(lo.x+j0+ix, lo.y+l +iz, 0), + &jx_arr(lo.x+j0+ix, lo.y+l +iz, 0), sx0[ix]*sz [iz]*wqx); amrex::Gpu::Atomic::Add( - &jy_arr(lo.x+j +ix, lo.y+l +iz, 0), + &jy_arr(lo.x+j +ix, lo.y+l +iz, 0), sx [ix]*sz [iz]*wqy); amrex::Gpu::Atomic::Add( - &jz_arr(lo.x+j +ix, lo.y+l0+iz, 0), + &jz_arr(lo.x+j +ix, lo.y+l0+iz, 0), sx [ix]*sz0[iz]*wqz); } } @@ -156,7 +156,7 @@ void doDepositionShapeN(const amrex::Real * const xp, &jx_arr(lo.x+j0+ix, lo.y+k +iy, lo.z+l +iz), sx0[ix]*sy [iy]*sz [iz]*wqx); amrex::Gpu::Atomic::Add( - &jy_arr(lo.x+j +ix, lo.y+k0+iy, lo.z+l +iz), + &jy_arr(lo.x+j +ix, lo.y+k0+iy, lo.z+l +iz), sx [ix]*sy0[iy]*sz [iz]*wqy); amrex::Gpu::Atomic::Add( &jz_arr(lo.x+j +ix, lo.y+k +iy, lo.z+l0+iz), @@ -240,7 +240,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c; // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr - amrex::ParallelFor( + amrex::ParallelFor( np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { @@ -302,7 +302,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, const amrex::Real x_new = (xp[ip] - xmin)*dxi; const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif -#if (defined WARPX_DIM_3D) +#if (defined WARPX_DIM_3D) const amrex::Real y_new = (yp[ip] - ymin)*dyi; const amrex::Real y_old = y_new - dtsdy0*uyp[ip]*gaminv; #endif @@ -336,7 +336,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, #if (defined WARPX_DIM_3D) const int j_new = compute_shape_factor<depos_order>(sy_new+1, y_new); const int j_old = compute_shifted_shape_factor<depos_order>(sy_old, y_old, j_new); -#endif +#endif const int k_new = compute_shape_factor<depos_order>(sz_new+1, z_new); const int k_old = compute_shifted_shape_factor<depos_order>(sz_old, z_old, k_new); @@ -418,7 +418,7 @@ void doEsirkepovDepositionShapeN (const amrex::Real * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -2.*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode* + const Complex djt_cmplx = -2.*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(double)imode* (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index d8d1d78ef..6727b0aa9 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -97,7 +97,7 @@ void doGatherShapeN(const amrex::Real * const xp, Bxp[ip] = 0; Byp[ip] = 0; Bzp[ip] = 0; - // Each field is gathered in a separate block of + // Each field is gathered in a separate block of // AMREX_SPACEDIM nested loops because the deposition // order can differ for each component of each field // when lower_in_v is set to 1 diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 305baf4ff..ac261b177 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -51,24 +51,24 @@ public: const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks); /// - /// This evolves all the particles by one PIC time step, including charge deposition, the + /// This evolves all the particles by one PIC time step, including charge deposition, the /// field solve, and pushing the particles, for all the species in the MultiParticleContainer. /// This is the electrostatic version. /// - void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, - amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, + void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, + amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, amrex::Real t, amrex::Real dt); /// - /// This pushes the particle positions by one half time step for all the species in the + /// This pushes the particle positions by one half time step for all the species in the /// MultiParticleContainer. It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. This is the electrostatic version. - /// + /// void PushXES (amrex::Real dt); /// /// This deposits the particle charge onto rho, accumulating the value for all the species - /// in the MultiParticleContainer. rho is assumed to contain node-centered multifabs. + /// in the MultiParticleContainer. rho is assumed to contain node-centered multifabs. /// This version is hard-coded for CIC deposition. /// void DepositCharge(amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, bool local = false); @@ -79,7 +79,7 @@ public: /// amrex::Real sumParticleCharge(bool local = false); #endif // WARPX_DO_ELECTROSTATIC - + /// /// Performs the field gather operation using the input fields E and B, for all the species /// in the MultiParticleContainer. This is the electromagnetic version of the field gather. @@ -87,7 +87,7 @@ public: void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, - const amrex::MultiFab& By, const amrex::MultiFab& Bz); + const amrex::MultiFab& By, const amrex::MultiFab& Bz); /// /// This evolves all the particles by one PIC time step, including current deposition, the @@ -97,8 +97,8 @@ public: void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, - amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, - amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, + amrex::MultiFab& jx, amrex::MultiFab& jy, amrex::MultiFab& jz, + amrex::MultiFab* cjx, amrex::MultiFab* cjy, amrex::MultiFab* cjz, amrex::MultiFab* rho, amrex::MultiFab* crho, const amrex::MultiFab* cEx, const amrex::MultiFab* cEy, const amrex::MultiFab* cEz, const amrex::MultiFab* cBx, const amrex::MultiFab* cBy, const amrex::MultiFab* cBz, @@ -108,7 +108,7 @@ public: /// This pushes the particle positions by one half time step for all the species in the /// MultiParticleContainer. It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. This is the electromagnetic version. - /// + /// void PushX (amrex::Real dt); /// @@ -116,7 +116,7 @@ public: /// MultiParticleContainer. It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. It is also used to synchronize particles at the /// the end of the run. This is the electromagnetic version. - /// + /// void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz); @@ -124,7 +124,7 @@ public: /// /// This deposits the particle charge onto a node-centered MultiFab and returns a unique ptr /// to it. The charge density is accumulated over all the particles in the MultiParticleContainer - /// + /// std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false); void doFieldIonization (); @@ -132,7 +132,7 @@ public: void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; - + void Restart (const std::string& dir); void PostRestart (); @@ -178,7 +178,7 @@ public: const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, amrex::Vector<WarpXParticleContainer::DiagnosticParticleData>& parts) const; - // Inject particles during the simulation (for particles entering the + // Inject particles during the simulation (for particles entering the // simulation domain after some iterations, due to flowing plasma and/or // moving window). void ContinuousInjection(const amrex::RealBox& injection_box) const; @@ -189,7 +189,7 @@ public: std::vector<std::string> GetSpeciesNames() const { return species_names; } PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } - + protected: // Particle container types @@ -218,10 +218,10 @@ private: void mapSpeciesProduct (); int getSpeciesID (std::string product_str); - + // Number of species dumped in BoostedFrameDiagnostics int nspecies_boosted_frame_diags = 0; - // map_species_boosted_frame_diags[i] is the species ID in + // map_species_boosted_frame_diags[i] is the species ID in // MultiParticleContainer for 0<i<nspecies_boosted_frame_diags std::vector<int> map_species_boosted_frame_diags; int do_boosted_frame_diags = 0; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 640439337..612583e2f 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -23,7 +23,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[i]->m_deposit_on_main_grid = m_deposit_on_main_grid[i]; allcontainers[i]->m_gather_from_main_grid = m_gather_from_main_grid[i]; } - + for (int i = nspecies; i < nspecies+nlasers; ++i) { allcontainers[i].reset(new LaserParticleContainer(amr_core, i, lasers_names[i-nspecies])); } @@ -359,7 +359,7 @@ MultiParticleContainer { BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - + // Loop over particle species for (int i = 0; i < nspecies_boosted_frame_diags; ++i){ int isp = map_species_boosted_frame_diags[i]; @@ -368,41 +368,41 @@ MultiParticleContainer pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. - + // Loop over AMR levels for (int lev = 0; lev <= pc->finestLevel(); ++lev){ - // Loop over [grid index][tile index] pairs - // and Fills parts[species number i] with particle data from all grids and - // tiles in diagnostic_particles. parts contains particles from all + // Loop over [grid index][tile index] pairs + // and Fills parts[species number i] with particle data from all grids and + // tiles in diagnostic_particles. parts contains particles from all // AMR levels indistinctly. for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ // it->first is the [grid index][tile index] key - // it->second is the corresponding + // it->second is the corresponding // WarpXParticleContainer::DiagnosticParticleData value parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); - + parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), it->second.GetRealData(DiagIdx::x ).begin(), it->second.GetRealData(DiagIdx::x ).end()); - + parts[i].GetRealData(DiagIdx::y).insert( parts[i].GetRealData(DiagIdx::y ).end(), it->second.GetRealData(DiagIdx::y ).begin(), it->second.GetRealData(DiagIdx::y ).end()); - + parts[i].GetRealData(DiagIdx::z).insert( parts[i].GetRealData(DiagIdx::z ).end(), it->second.GetRealData(DiagIdx::z ).begin(), it->second.GetRealData(DiagIdx::z ).end()); - + parts[i].GetRealData(DiagIdx::ux).insert( parts[i].GetRealData(DiagIdx::ux).end(), it->second.GetRealData(DiagIdx::ux).begin(), it->second.GetRealData(DiagIdx::ux).end()); - + parts[i].GetRealData(DiagIdx::uy).insert( parts[i].GetRealData(DiagIdx::uy).end(), it->second.GetRealData(DiagIdx::uy).begin(), it->second.GetRealData(DiagIdx::uy).end()); - + parts[i].GetRealData(DiagIdx::uz).insert( parts[i].GetRealData(DiagIdx::uz).end(), it->second.GetRealData(DiagIdx::uz).begin(), it->second.GetRealData(DiagIdx::uz).end()); @@ -413,7 +413,7 @@ MultiParticleContainer /* \brief Continuous injection for particles initially outside of the domain. * \param injection_box: Domain where new particles should be injected. - * Loop over all WarpXParticleContainer in MultiParticleContainer and + * Loop over all WarpXParticleContainer in MultiParticleContainer and * calls virtual function ContinuousInjection. */ void @@ -429,7 +429,7 @@ MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const /* \brief Update position of continuous injection parameters. * \param dt: simulation time step (level 0) - * All classes inherited from WarpXParticleContainer do not have + * All classes inherited from WarpXParticleContainer do not have * a position to update (PhysicalParticleContainer does not do anything). */ void @@ -457,7 +457,7 @@ MultiParticleContainer::doContinuousInjection () const } /* \brief Get ID of product species of each species. - * The users specifies the name of the product species, + * The users specifies the name of the product species, * this routine get its ID. */ void @@ -465,8 +465,8 @@ MultiParticleContainer::mapSpeciesProduct () { for (int i=0; i<nspecies; i++){ auto& pc = allcontainers[i]; - // If species pc has ionization on, find species with name - // pc->ionization_product_name and store its ID into + // If species pc has ionization on, find species with name + // pc->ionization_product_name and store its ID into // pc->ionization_product. if (pc->do_field_ionization){ int i_product = getSpeciesID(pc->ionization_product_name); @@ -577,7 +577,7 @@ namespace } // --- product runtime attribs GpuArray<Real*,6> runtime_attribs_product; - bool do_boosted_product = WarpX::do_boosted_frame_diagnostic + bool do_boosted_product = WarpX::do_boosted_frame_diagnostic && pc_product->DoBoostedFrameDiags(); if (do_boosted_product) { std::map<std::string, int> comps_product = pc_product->getParticleComps(); @@ -624,7 +624,7 @@ namespace attribs_product[ia][ip] = attribs_source[ia][is]; } // Update xold etc. if boosted frame diagnostics required - // for product species. Fill runtime attribs with a copy of + // for product species. Fill runtime attribs with a copy of // current properties (xold = x etc.). if (do_boosted_product) { runtime_attribs_product[0][ip] = p_source.pos(0); @@ -647,7 +647,7 @@ MultiParticleContainer::doFieldIonization () // Loop over all species. // Ionized particles in pc_source create particles in pc_product for (auto& pc_source : allcontainers){ - + // Skip if not ionizable if (!pc_source->do_field_ionization){ continue; } @@ -661,7 +661,7 @@ MultiParticleContainer::doFieldIonization () // they do not exist (or if they were defined by default, i.e., // without runtime component). #ifdef _OPENMP - // Touch all tiles of source species in serial if runtime attribs + // Touch all tiles of source species in serial if runtime attribs for (MFIter mfi = pc_source->MakeMFIter(lev); mfi.isValid(); ++mfi) { const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); @@ -683,7 +683,7 @@ MultiParticleContainer::doFieldIonization () MFItInfo info; if (pc_source->do_tiling && Gpu::notInLaunchRegion()) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - pc_product->do_tiling, + pc_product->do_tiling, "For ionization, either all or none of the " "particle species must use tiling."); info.EnableTiling(pc_source->tile_size); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c7494fbdf..7946b4650 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -24,7 +24,7 @@ public: void InitIonizationModule (); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks) override; @@ -32,7 +32,7 @@ public: amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC - + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -99,14 +99,14 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} void SplitParticles(int lev); - + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) override; @@ -125,7 +125,7 @@ public: std::array<amrex::Real, 3> u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 318ad4664..99c6973f9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -51,7 +51,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); @@ -68,9 +68,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -142,7 +142,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -154,7 +154,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution<double> distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -387,11 +387,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector<int> cellid_v; if (refine_injection and lev == 0) @@ -446,7 +446,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray<Real,AMREX_SPACEDIM> overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -458,9 +458,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -627,7 +627,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4<Real> const& costarr = cost->array(mfi); @@ -865,7 +865,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -915,7 +915,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -949,11 +949,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // Get instances of NCI Godfrey filters const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; @@ -978,9 +978,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1183,7 +1183,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1199,14 +1199,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1218,7 +1218,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); @@ -1234,7 +1234,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1272,13 +1272,13 @@ PhysicalParticleContainer::Evolve (int lev, eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); ceyfab = &filtered_Ey; - + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); cbxfab = &filtered_Bx; - + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); @@ -1286,14 +1286,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1303,7 +1303,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1317,7 +1317,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1337,7 +1337,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1406,7 +1406,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; i<np; i++){ auto& p = particles[i]; if (p.id() == DoSplitParticleID){ - // If particle is tagged, split it and put the + // If particle is tagged, split it and put the // split particles in local arrays psplit_x etc. np_split_to_add += np_split; #if (AMREX_SPACEDIM==2) @@ -1503,11 +1503,11 @@ PhysicalParticleContainer::SplitParticles(int lev) } // Add local arrays psplit_x etc. to the temporary // particle container pctmp_split. Split particles - // are tagged with p.id()=NoSplitParticleID so that + // are tagged with p.id()=NoSplitParticleID so that // they are not re-split when entering a higher level // AddNParticles calls Redistribute, so that particles // in pctmp_split are in the proper grids and tiles - pctmp_split.AddNParticles(lev, + pctmp_split.AddNParticles(lev, np_split_to_add, psplit_x.dataPtr(), psplit_y.dataPtr(), @@ -1559,12 +1559,12 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, if (do_field_ionization){ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr(); } - + // Loop over the particles and update their momentum const Real q = this->charge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1613,7 +1613,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1653,7 +1653,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -1701,7 +1701,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1717,7 +1717,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1757,7 +1757,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1804,7 +1804,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1839,7 +1839,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1864,7 +1864,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1874,7 +1874,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1901,14 +1901,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1918,26 +1918,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4<const Real>& ex_arr = exfab->array(); const Array4<const Real>& ey_arr = eyfab->array(); const Array4<const Real>& ez_arr = ezfab->array(); const Array4<const Real>& bx_arr = bxfab->array(); const Array4<const Real>& by_arr = byfab->array(); const Array4<const Real>& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2020,10 +2020,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2038,18 +2038,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2095,7 +2095,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index b920ece0a..a635e9aa0 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -60,7 +60,7 @@ public: virtual void ReadHeader (std::istream& is) override; virtual void WriteHeader (std::ostream& os) const override; - + private: // User input quantities diff --git a/Source/Particles/ShapeFactors.H b/Source/Particles/ShapeFactors.H index 9d185714a..be79a4871 100644 --- a/Source/Particles/ShapeFactors.H +++ b/Source/Particles/ShapeFactors.H @@ -1,7 +1,7 @@ #ifndef SHAPEFACTORS_H_ #define SHAPEFACTORS_H_ -// Compute shape factor and return index of leftmost cell where +// Compute shape factor and return index of leftmost cell where // particle writes. // Specialized templates are defined below for orders 0 to 3. template <int depos_order> diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index f7e46b2d3..06540871b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -73,19 +73,19 @@ public: const amrex::Cuda::ManagedDeviceVector<amrex::Real>& y, const amrex::Cuda::ManagedDeviceVector<amrex::Real>& z); #endif - const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { - return GetStructOfArrays().GetRealData(); + const std::array<RealVector, PIdx::nattribs>& GetAttribs () const { + return GetStructOfArrays().GetRealData(); } - - std::array<RealVector, PIdx::nattribs>& GetAttribs () { - return GetStructOfArrays().GetRealData(); + + std::array<RealVector, PIdx::nattribs>& GetAttribs () { + return GetStructOfArrays().GetRealData(); } - const RealVector& GetAttribs (int comp) const { + const RealVector& GetAttribs (int comp) const { return GetStructOfArrays().GetRealData(comp); } - - RealVector& GetAttribs (int comp) { + + RealVector& GetAttribs (int comp) { return GetStructOfArrays().GetRealData(comp); } @@ -102,15 +102,15 @@ class WarpXParticleContainer public: friend MultiParticleContainer; - // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components // and 0 int components for the particle data. using DiagnosticParticleData = amrex::StructOfArrays<DiagIdx::nattribs, 0>; // DiagnosticParticles is a vector, with one element per MR level. - // DiagnosticParticles[lev] is typically a key-value pair where the key is - // a pair [grid_index, tile_index], and the value is the corresponding + // DiagnosticParticles[lev] is typically a key-value pair where the key is + // a pair [grid_index, tile_index], and the value is the corresponding // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector<std::map<std::pair<int, int>, DiagnosticParticleData> >; - + WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); virtual ~WarpXParticleContainer() {} @@ -124,12 +124,12 @@ public: const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) {} -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void EvolveES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, 3> >& E, - amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, + amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, amrex::Real t, amrex::Real dt) = 0; #endif // WARPX_DO_ELECTROSTATIC - + virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, @@ -143,10 +143,10 @@ public: virtual void PostRestart () = 0; virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) {} - + void AllocData (); /// @@ -154,7 +154,7 @@ public: /// It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. /// This is the electrostatic version of the particle push. - /// + /// void PushXES (amrex::Real dt); /// @@ -162,13 +162,13 @@ public: /// It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. /// This is the electromagnetic version of the particle push. - /// + /// void PushX ( amrex::Real dt); void PushX (int lev, amrex::Real dt); /// /// This pushes the particle momenta by dt. - /// + /// virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -208,9 +208,9 @@ public: int depos_lev, amrex::Real dt); - // If particles start outside of the domain, ContinuousInjection - // makes sure that they are initialized when they enter the domain, and - // NOT before. Virtual function, overriden by derived classes. + // If particles start outside of the domain, ContinuousInjection + // makes sure that they are initialized when they enter the domain, and + // NOT before. Virtual function, overriden by derived classes. // Current status: // PhysicalParticleContainer: implemented. // LaserParticleContainer: implemented. @@ -219,7 +219,7 @@ public: // Update optional sub-class-specific injection location. virtual void UpdateContinuousInjectionPosition(amrex::Real dt) {} - /// + /// /// This returns the total charge for all the particles in this ParticleContainer. /// This is needed when solving Poisson's equation with periodic boundary conditions. /// @@ -259,9 +259,9 @@ public: // split along axes (0) or diagonals (1) int split_type = 0; - using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; + using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddIntComp; - + void AddRealComp (const std::string& name, bool comm=true) { particle_comps[name] = NumRealComps(); @@ -274,7 +274,7 @@ public: AddIntComp(comm); } - int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) @@ -286,7 +286,7 @@ protected: std::map<std::string, int> particle_comps; std::map<std::string, int> particle_icomps; - + int species_id; amrex::Real charge; @@ -300,9 +300,9 @@ protected: static int do_not_push; - // Whether to allow particles outside of the simulation domain to be + // Whether to allow particles outside of the simulation domain to be // initialized when they enter the domain. - // This is currently required because continuous injection does not + // This is currently required because continuous injection does not // support all features allowed by direct injection. int do_continuous_injection = 0; @@ -316,7 +316,7 @@ protected: amrex::Gpu::ManagedVector<amrex::Real> adk_prefactor; amrex::Gpu::ManagedVector<amrex::Real> adk_exp_prefactor; std::string physical_element; - + int do_boosted_frame_diags = 1; amrex::Vector<amrex::FArrayBox> local_rho; @@ -326,20 +326,20 @@ protected: using DataContainer = amrex::Gpu::ManagedDeviceVector<amrex::Real>; using PairIndex = std::pair<int, int>; - + amrex::Vector<DataContainer> m_xp, m_yp, m_zp, m_giv; - // Whether to dump particle quantities. + // Whether to dump particle quantities. // If true, particle position is always dumped. int plot_species = 1; - // For all particle attribs (execept position), whether or not + // For all particle attribs (execept position), whether or not // to dump to file. amrex::Vector<int> plot_flags; // list of names of attributes to dump. amrex::Vector<std::string> plot_vars; - + amrex::Vector<std::map<PairIndex, std::array<DataContainer, TmpIdx::nattribs> > > tmp_particle_data; - + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 4e80374d8..83de9fb0e 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -228,7 +228,7 @@ WarpXParticleContainer::AddNParticles (int lev, #endif if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); } particle_tile.push_back(p); @@ -242,7 +242,7 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); if ( (NumRuntimeRealComps()>0) || (NumRuntimeIntComps()>0) ){ - auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) @@ -289,7 +289,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, Real dt) -{ +{ AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || (depos_lev==(lev )), "Deposition buffers only work for lev-1"); @@ -316,7 +316,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(depos_lev); tilebox = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Staggered tile boxes (different in each direction) Box tbx = convert(tilebox, WarpX::jx_nodal_flag); Box tby = convert(tilebox, WarpX::jy_nodal_flag); @@ -360,11 +360,11 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, // Lower corner of tile box physical domain // Note that this includes guard cells since it is after tilebox.ngrow const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(tilebox, depos_lev); - // xyzmin is built on pti.tilebox(), so it does + // xyzmin is built on pti.tilebox(), so it does // not include staggering, so the stagger_shift has to be done by hand. - // Alternatively, we could define xyzminx from tbx (and the same for 3 + // Alternatively, we could define xyzminx from tbx (and the same for 3 // directions and for jx, jy, jz). This way, sx0 would not be needed. - // Better for memory? worth trying? + // Better for memory? worth trying? const Dim3 lo = lbound(tilebox); BL_PROFILE_VAR_START(blp_deposit); @@ -431,7 +431,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, ion_lev is a null pointer. * \param rho : Full array of charge density * \param icomp : Component of rho into which charge is deposited. - 0: old value (before particle push). + 0: old value (before particle push). 1: new value (after particle push). * \param offset : Index of first particle for which charge is deposited * \param np_to_depose: Number of particles for which charge is deposited. @@ -471,7 +471,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const IntVect& ref_ratio = WarpX::RefRatio(depos_lev); tilebox = amrex::coarsen(pti.tilebox(),ref_ratio); } - + tilebox.grow(ngRho); const int nc = (rho->nComp() == 1 ? 1 : rho->nComp()/2); diff --git a/Source/Particles/deposit_cic.F90 b/Source/Particles/deposit_cic.F90 index e4e1ace03..1fe80016f 100644 --- a/Source/Particles/deposit_cic.F90 +++ b/Source/Particles/deposit_cic.F90 @@ -7,11 +7,11 @@ module warpx_ES_deposit_cic contains -! This routine computes the charge density due to the particles using cloud-in-cell +! This routine computes the charge density due to the particles using cloud-in-cell ! deposition. The Fab rho is assumed to be node-centered. ! ! Arguments: -! particles : a pointer to the particle array-of-structs +! particles : a pointer to the particle array-of-structs ! ns : the stride length of particle struct (the size of the struct in number of reals) ! np : the number of particles ! weights : the particle weights diff --git a/Source/Particles/interpolate_cic.F90 b/Source/Particles/interpolate_cic.F90 index bc2c4f37e..005ab35f4 100644 --- a/Source/Particles/interpolate_cic.F90 +++ b/Source/Particles/interpolate_cic.F90 @@ -12,7 +12,7 @@ contains ! node-centered. ! ! Arguments: -! particles : a pointer to the particle array-of-structs +! particles : a pointer to the particle array-of-structs ! ns : the stride length of particle struct (the size of the struct in number of reals) ! np : the number of particles ! Ex_p : the electric field in the x-direction at the particle positions (output) @@ -91,7 +91,7 @@ contains wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) - + end do end subroutine warpx_interpolate_cic_3d @@ -184,7 +184,7 @@ contains lx = (particles(1, n) - plo(1))*inv_dx(1) ly = (particles(2, n) - plo(2))*inv_dx(2) lz = (particles(3, n) - plo(3))*inv_dx(3) - + i = floor(lx) j = floor(ly) k = floor(lz) @@ -203,11 +203,11 @@ contains wx_hi = lx - i wy_hi = ly - j wz_hi = lz - k - + wx_lo = 1.0d0 - wx_hi wy_lo = 1.0d0 - wy_hi wz_lo = 1.0d0 - wz_hi - + Ex_p(n) = wx_lo*wy_lo*wz_lo*cEx(i, j, k ) + & wx_lo*wy_lo*wz_hi*cEx(i, j, k+1) + & wx_lo*wy_hi*wz_lo*cEx(i, j+1, k ) + & @@ -225,7 +225,7 @@ contains wx_hi*wy_lo*wz_hi*cEy(i+1, j, k+1) + & wx_hi*wy_hi*wz_lo*cEy(i+1, j+1, k ) + & wx_hi*wy_hi*wz_hi*cEy(i+1, j+1, k+1) - + Ez_p(n) = wx_lo*wy_lo*wz_lo*cEz(i, j, k ) + & wx_lo*wy_lo*wz_hi*cEz(i, j, k+1) + & wx_lo*wy_hi*wz_lo*cEz(i, j+1, k ) + & @@ -234,14 +234,14 @@ contains wx_hi*wy_lo*wz_hi*cEz(i+1, j, k+1) + & wx_hi*wy_hi*wz_lo*cEz(i+1, j+1, k ) + & wx_hi*wy_hi*wz_hi*cEz(i+1, j+1, k+1) - + ! otherwise use the fine else wx_hi = lx - i wy_hi = ly - j wz_hi = lz - k - + wx_lo = 1.0d0 - wx_hi wy_lo = 1.0d0 - wy_hi wz_lo = 1.0d0 - wz_hi @@ -254,7 +254,7 @@ contains wx_hi*wy_lo*wz_hi*Ex(i+1, j, k+1) + & wx_hi*wy_hi*wz_lo*Ex(i+1, j+1, k ) + & wx_hi*wy_hi*wz_hi*Ex(i+1, j+1, k+1) - + Ey_p(n) = wx_lo*wy_lo*wz_lo*Ey(i, j, k ) + & wx_lo*wy_lo*wz_hi*Ey(i, j, k+1) + & wx_lo*wy_hi*wz_lo*Ey(i, j+1, k ) + & @@ -263,7 +263,7 @@ contains wx_hi*wy_lo*wz_hi*Ey(i+1, j, k+1) + & wx_hi*wy_hi*wz_lo*Ey(i+1, j+1, k ) + & wx_hi*wy_hi*wz_hi*Ey(i+1, j+1, k+1) - + Ez_p(n) = wx_lo*wy_lo*wz_lo*Ez(i, j, k ) + & wx_lo*wy_lo*wz_hi*Ez(i, j, k+1) + & wx_lo*wy_hi*wz_lo*Ez(i, j+1, k ) + & @@ -272,7 +272,7 @@ contains wx_hi*wy_lo*wz_hi*Ez(i+1, j, k+1) + & wx_hi*wy_hi*wz_lo*Ez(i+1, j+1, k ) + & wx_hi*wy_hi*wz_hi*Ez(i+1, j+1, k+1) - + end if end do @@ -314,7 +314,7 @@ contains lx = (particles(1, n) - plo(1))*inv_dx(1) ly = (particles(2, n) - plo(2))*inv_dx(2) - + i = floor(lx) j = floor(ly) @@ -329,10 +329,10 @@ contains wx_hi = lx - i wy_hi = ly - j - + wx_lo = 1.0d0 - wx_hi wy_lo = 1.0d0 - wy_hi - + Ex_p(n) = wx_lo*wy_lo*cEx(i, j ) + & wx_lo*wy_hi*cEx(i, j+1) + & wx_hi*wy_lo*cEx(i+1, j ) + & @@ -342,13 +342,13 @@ contains wx_lo*wy_hi*cEy(i, j+1) + & wx_hi*wy_lo*cEy(i+1, j ) + & wx_hi*wy_hi*cEy(i+1, j+1) - + ! otherwise use the fine else wx_hi = lx - i wy_hi = ly - j - + wx_lo = 1.0d0 - wx_hi wy_lo = 1.0d0 - wy_hi @@ -356,7 +356,7 @@ contains wx_lo*wy_hi*Ex(i, j+1) + & wx_hi*wy_lo*Ex(i+1, j ) + & wx_hi*wy_hi*Ex(i+1, j+1) - + Ey_p(n) = wx_lo*wy_lo*Ey(i, j ) + & wx_lo*wy_hi*Ey(i, j+1) + & wx_hi*wy_lo*Ey(i+1, j ) + & diff --git a/Source/Particles/push_particles_ES.F90 b/Source/Particles/push_particles_ES.F90 index 9c7bd9e92..53dd3c181 100644 --- a/Source/Particles/push_particles_ES.F90 +++ b/Source/Particles/push_particles_ES.F90 @@ -8,13 +8,13 @@ module warpx_ES_push_particles contains ! -! This routine updates the particle positions and velocities using the +! This routine updates the particle positions and velocities using the ! leapfrog time integration algorithm, given the electric fields at the ! particle positions. It also enforces specular reflection off the domain ! walls. ! ! Arguments: -! particles : a pointer to the particle array-of-structs +! particles : a pointer to the particle array-of-structs ! ns : the stride length of particle struct (the size of the struct in number of reals) ! np : the number of particles ! vx_p : the particle x-velocities @@ -30,7 +30,7 @@ contains ! prob_hi : the right-hand corner of the problem domain ! subroutine warpx_push_leapfrog_3d(particles, ns, np, & - vx_p, vy_p, vz_p, & + vx_p, vy_p, vz_p, & Ex_p, Ey_p, Ez_p, & charge, mass, dt, & prob_lo, prob_hi) & @@ -43,7 +43,7 @@ contains real(amrex_real), intent(in) :: mass real(amrex_real), intent(in) :: dt real(amrex_real), intent(in) :: prob_lo(3), prob_hi(3) - + integer n real(amrex_real) fac @@ -69,7 +69,7 @@ contains vx_p(n) = -vx_p(n) end do -! ... y... +! ... y... do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) if (particles(2, n) .lt. prob_lo(2)) then particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) @@ -107,7 +107,7 @@ contains real(amrex_real), intent(in) :: mass real(amrex_real), intent(in) :: dt real(amrex_real), intent(in) :: prob_lo(2), prob_hi(2) - + integer n real(amrex_real) fac @@ -131,7 +131,7 @@ contains vx_p(n) = -vx_p(n) end do -! ... y... +! ... y... do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) if (particles(2, n) .lt. prob_lo(2)) then particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) @@ -152,7 +152,7 @@ contains ! from the velocities after particle initialization. ! ! Arguments: -! particles : a pointer to the particle array-of-structs +! particles : a pointer to the particle array-of-structs ! ns : the stride length of particle struct (the size of the struct in number of reals) ! np : the number of particles ! xx_p : the electric field in the x-direction at the particle positions @@ -190,7 +190,7 @@ contains vx_p(n) = -vx_p(n) end do -! ... y... +! ... y... do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) if (particles(2, n) .lt. prob_lo(2)) then particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) @@ -209,7 +209,7 @@ contains end if vz_p(n) = -vz_p(n) end do - + end do end subroutine warpx_push_leapfrog_positions_3d @@ -241,7 +241,7 @@ contains vx_p(n) = -vx_p(n) end do -! ... y... +! ... y... do while (particles(2, n) .lt. prob_lo(2) .or. particles(2, n) .gt. prob_hi(2)) if (particles(2, n) .lt. prob_lo(2)) then particles(2, n) = 2.d0*prob_lo(2) - particles(2, n) diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 2fb7200bb..b59f80495 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -8,7 +8,7 @@ #include <WarpXUtil.H> #include <WarpX_py.H> -namespace +namespace { double** getMultiFabPointers(const amrex::MultiFab& mf, int *num_boxes, int *ncomps, int *ngrow, int **shapes) { @@ -19,7 +19,7 @@ namespace if (mf.nComp() > 1) shapesize += 1; *shapes = (int*) malloc(shapesize * (*num_boxes) * sizeof(int)); double** data = (double**) malloc((*num_boxes) * sizeof(double*)); - + #ifdef _OPENMP #pragma omp parallel #endif @@ -27,7 +27,7 @@ namespace int i = mfi.LocalIndex(); data[i] = (double*) mf[mfi].dataPtr(); for (int j = 0; j < AMREX_SPACEDIM; ++j) { - (*shapes)[shapesize*i+j] = mf[mfi].box().length(j); + (*shapes)[shapesize*i+j] = mf[mfi].box().length(j); } if (mf.nComp() > 1) (*shapes)[shapesize*i+AMREX_SPACEDIM] = mf.nComp(); } @@ -69,12 +69,12 @@ extern "C" return WarpX::l_lower_order_in_v; } - int warpx_nComps() + int warpx_nComps() { - return PIdx::nattribs; + return PIdx::nattribs; } - int warpx_SpaceDim() + int warpx_SpaceDim() { return AMREX_SPACEDIM; } @@ -87,7 +87,7 @@ extern "C" #ifdef BL_USE_MPI void amrex_init_with_inited_mpi (int argc, char* argv[], MPI_Comm mpicomm) { - amrex::Initialize(argc,argv,true,mpicomm); + amrex::Initialize(argc,argv,true,mpicomm); } #endif diff --git a/Source/Python/WarpXWrappers.h b/Source/Python/WarpXWrappers.h index 0a476b5fc..233a0b930 100644 --- a/Source/Python/WarpXWrappers.h +++ b/Source/Python/WarpXWrappers.h @@ -26,9 +26,9 @@ extern "C" { #endif void amrex_finalize (int finalize_mpi); - + void warpx_init (); - + void warpx_finalize (); typedef void(*WARPX_CALLBACK_PY_FUNC_0)(); @@ -47,41 +47,41 @@ extern "C" { void warpx_set_callback_py_appliedfields (WARPX_CALLBACK_PY_FUNC_0); void warpx_evolve (int numsteps); // -1 means the inputs parameter will be used. - + void warpx_addNParticles(int speciesnumber, int lenx, double* x, double* y, double* z, double* vx, double* vy, double* vz, int nattr, double* attr, int uniqueparticles); void warpx_ConvertLabParamsToBoost(); - + double warpx_getProbLo(int dir); - + double warpx_getProbHi(int dir); - + long warpx_getNumParticles(int speciesnumber); - - double** warpx_getEfield(int lev, int direction, + + double** warpx_getEfield(int lev, int direction, int *return_size, int* ncomps, int* ngrow, int **shapes); - - int* warpx_getEfieldLoVects(int lev, int direction, + + int* warpx_getEfieldLoVects(int lev, int direction, int *return_size, int* ngrow); - - double** warpx_getBfield(int lev, int direction, + + double** warpx_getBfield(int lev, int direction, int *return_size, int* ncomps, int* ngrow, int **shapes); - - int* warpx_getBfieldLoVects(int lev, int direction, + + int* warpx_getBfieldLoVects(int lev, int direction, int *return_size, int* ngrow); - - double** warpx_getCurrentDensity(int lev, int direction, + + double** warpx_getCurrentDensity(int lev, int direction, int *return_size, int* ncomps, int* ngrow, int **shapes); - - int* warpx_getCurrentDensityLoVects(int lev, int direction, + + int* warpx_getCurrentDensityLoVects(int lev, int direction, int *return_size, int* ngrow); - + double** warpx_getParticleStructs(int speciesnumber, int lev, int* num_tiles, int** particles_per_tile); - + double** warpx_getParticleArrays(int speciesnumber, int comp, int lev, int* num_tiles, int** particles_per_tile); diff --git a/Source/Utils/IonizationEnergiesTable.H b/Source/Utils/IonizationEnergiesTable.H index a7ab8210f..4c2898b91 100644 --- a/Source/Utils/IonizationEnergiesTable.H +++ b/Source/Utils/IonizationEnergiesTable.H @@ -32,105 +32,105 @@ std::map<std::string, int> ion_map_ids = { const int nelements = 22; const int ion_atomic_numbers[nelements] = { - 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, - 11, 12, 13, 14, 15, 16, 17, 18, 36, 37, + 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, + 11, 12, 13, 14, 15, 16, 17, 18, 36, 37, 54, 86}; const int ion_energy_offsets[nelements] = { - 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, - 55, 66, 78, 91, 105, 120, 136, 153, 171, 207, + 0, 1, 3, 6, 10, 15, 21, 28, 36, 45, + 55, 66, 78, 91, 105, 120, 136, 153, 171, 207, 244, 298}; const int energies_tab_length = 384; const amrex::Real table_ionization_energies[energies_tab_length]{ // H - 13.59843449, + 13.59843449, // He - 24.58738880, 54.4177650, + 24.58738880, 54.4177650, // Li - 5.39171495, 75.6400964, 122.4543581, + 5.39171495, 75.6400964, 122.4543581, // Be - 9.322699, 18.21115, 153.896203, 217.7185843, + 9.322699, 18.21115, 153.896203, 217.7185843, // B - 8.298019, 25.15483, 37.93058, 259.3715, 340.226020, + 8.298019, 25.15483, 37.93058, 259.3715, 340.226020, // C - 11.2602880, 24.383154, 47.88778, 64.49352, 392.090515, 489.993194, + 11.2602880, 24.383154, 47.88778, 64.49352, 392.090515, 489.993194, // N - 14.53413, 29.60125, 47.4453, 77.4735, 97.8901, 552.06732, 667.046116, + 14.53413, 29.60125, 47.4453, 77.4735, 97.8901, 552.06732, 667.046116, // O - 13.618055, 35.12112, 54.93554, 77.41350, 113.8990, 138.1189, 739.32682, - 871.40988, + 13.618055, 35.12112, 54.93554, 77.41350, 113.8990, 138.1189, 739.32682, + 871.40988, // F - 17.42282, 34.97081, 62.70798, 87.175, 114.249, 157.16311, 185.1868, - 953.89804, 1103.11747, + 17.42282, 34.97081, 62.70798, 87.175, 114.249, 157.16311, 185.1868, + 953.89804, 1103.11747, // Ne - 21.564540, 40.96297, 63.4233, 97.1900, 126.247, 157.934, 207.271, - 239.0970, 1195.80783, 1362.19915, + 21.564540, 40.96297, 63.4233, 97.1900, 126.247, 157.934, 207.271, + 239.0970, 1195.80783, 1362.19915, // Na - 5.1390769, 47.28636, 71.6200, 98.936, 138.404, 172.23, 208.504, - 264.192, 299.856, 1465.13449, 1648.70218, + 5.1390769, 47.28636, 71.6200, 98.936, 138.404, 172.23, 208.504, + 264.192, 299.856, 1465.13449, 1648.70218, // Mg - 7.646236, 15.035271, 80.1436, 109.2654, 141.33, 186.76, 225.02, - 265.924, 327.99, 367.489, 1761.80487, 1962.66365, + 7.646236, 15.035271, 80.1436, 109.2654, 141.33, 186.76, 225.02, + 265.924, 327.99, 367.489, 1761.80487, 1962.66365, // Al - 5.985769, 18.82855, 28.447642, 119.9924, 153.8252, 190.49, 241.76, - 284.64, 330.21, 398.65, 442.005, 2085.97700, 2304.14005, + 5.985769, 18.82855, 28.447642, 119.9924, 153.8252, 190.49, 241.76, + 284.64, 330.21, 398.65, 442.005, 2085.97700, 2304.14005, // Si - 8.15168, 16.34585, 33.49300, 45.14179, 166.767, 205.279, 246.57, - 303.59, 351.28, 401.38, 476.273, 523.415, 2437.65813, 2673.17753, + 8.15168, 16.34585, 33.49300, 45.14179, 166.767, 205.279, 246.57, + 303.59, 351.28, 401.38, 476.273, 523.415, 2437.65813, 2673.17753, // P - 10.486686, 19.76949, 30.20264, 51.44387, 65.02511, 220.430, 263.57, - 309.60, 372.31, 424.40, 479.44, 560.62, 611.741, 2816.90876, - 3069.8415, + 10.486686, 19.76949, 30.20264, 51.44387, 65.02511, 220.430, 263.57, + 309.60, 372.31, 424.40, 479.44, 560.62, 611.741, 2816.90876, + 3069.8415, // S - 10.36001, 23.33788, 34.86, 47.222, 72.5945, 88.0529, 280.954, - 328.794, 379.84, 447.7, 504.55, 564.41, 651.96, 706.994, - 3223.7807, 3494.1879, + 10.36001, 23.33788, 34.86, 47.222, 72.5945, 88.0529, 280.954, + 328.794, 379.84, 447.7, 504.55, 564.41, 651.96, 706.994, + 3223.7807, 3494.1879, // Cl - 12.967632, 23.81364, 39.80, 53.24, 67.68, 96.94, 114.2013, - 348.306, 400.851, 456.7, 530.0, 591.58, 656.30, 750.23, - 809.198, 3658.3437, 3946.2909, + 12.967632, 23.81364, 39.80, 53.24, 67.68, 96.94, 114.2013, + 348.306, 400.851, 456.7, 530.0, 591.58, 656.30, 750.23, + 809.198, 3658.3437, 3946.2909, // Ar - 15.7596117, 27.62967, 40.735, 59.58, 74.84, 91.290, 124.41, - 143.4567, 422.60, 479.76, 540.4, 619.0, 685.5, 755.13, - 855.5, 918.375, 4120.6656, 4426.2228, + 15.7596117, 27.62967, 40.735, 59.58, 74.84, 91.290, 124.41, + 143.4567, 422.60, 479.76, 540.4, 619.0, 685.5, 755.13, + 855.5, 918.375, 4120.6656, 4426.2228, // Kr - 13.9996053, 24.35984, 35.838, 50.85, 64.69, 78.49, 109.13, - 125.802, 233.0, 268, 308, 350, 391, 446, - 492, 540, 591, 640, 785, 831.6, 882.8, - 945, 999.0, 1042, 1155.0, 1205.23, 2928.9, 3072, - 3228, 3380, 3584, 3752.0, 3971, 4109.083, 17296.420, - 17936.209, + 13.9996053, 24.35984, 35.838, 50.85, 64.69, 78.49, 109.13, + 125.802, 233.0, 268, 308, 350, 391, 446, + 492, 540, 591, 640, 785, 831.6, 882.8, + 945, 999.0, 1042, 1155.0, 1205.23, 2928.9, 3072, + 3228, 3380, 3584, 3752.0, 3971, 4109.083, 17296.420, + 17936.209, // Rb - 4.1771280, 27.28954, 39.247, 52.20, 68.44, 82.9, 98.67, - 132.79, 150.628, 277.12, 313.1, 356.0, 400, 443, - 502, 550, 601, 654, 706.0, 857, 905.3, - 958.9, 1024, 1080, 1125, 1242.5, 1294.57, 3133.3, - 3281, 3443, 3600, 3815, 3988, 4214, 4356.865, - 18305.884, 18965.516, + 4.1771280, 27.28954, 39.247, 52.20, 68.44, 82.9, 98.67, + 132.79, 150.628, 277.12, 313.1, 356.0, 400, 443, + 502, 550, 601, 654, 706.0, 857, 905.3, + 958.9, 1024, 1080, 1125, 1242.5, 1294.57, 3133.3, + 3281, 3443, 3600, 3815, 3988, 4214, 4356.865, + 18305.884, 18965.516, // Xe - 12.1298436, 20.975, 31.05, 42.20, 54.1, 66.703, 91.6, - 105.9778, 179.84, 202.0, 229.02, 255.0, 281, 314, - 343, 374, 404, 434, 549, 582, 616, - 650, 700, 736, 818, 857.0, 1493, 1571, - 1653, 1742, 1826, 1919, 2023, 2113, 2209, - 2300, 2556, 2637, 2726, 2811, 2975, 3068, - 3243, 3333.8, 7660, 7889, 8144, 8382, 8971, - 9243, 9581, 9810.37, 40271.724, 41299.71, + 12.1298436, 20.975, 31.05, 42.20, 54.1, 66.703, 91.6, + 105.9778, 179.84, 202.0, 229.02, 255.0, 281, 314, + 343, 374, 404, 434, 549, 582, 616, + 650, 700, 736, 818, 857.0, 1493, 1571, + 1653, 1742, 1826, 1919, 2023, 2113, 2209, + 2300, 2556, 2637, 2726, 2811, 2975, 3068, + 3243, 3333.8, 7660, 7889, 8144, 8382, 8971, + 9243, 9581, 9810.37, 40271.724, 41299.71, // Rn - 10.74850, 21.4, 29.4, 36.9, 52.9, 64.0, 88.0, - 102.0, 154.0, 173.9, 195.0, 218.0, 240, 264, - 293, 317, 342, 367, 488, 520, 550, - 580, 640, 680, 760, 800, 850, 920, - 980, 1050, 1110, 1180, 1250, 1310, 1390, - 1460, 1520, 1590, 1660, 1720, 2033, 2094, - 2158, 2227, 2293, 2357, 2467, 2535, 2606, - 2674, 2944, 3010, 3082, 3149, 3433, 3510, - 3699, 3777, 6169, 6318, 6476, 6646, 6807, - 6964, 7283, 7450, 7630, 7800, 8260, 8410, - 8570, 8710, 9610, 9780, 10120, 10290, 21770, - 22160, 22600, 22990, 26310, 26830, 27490, 27903.1, + 10.74850, 21.4, 29.4, 36.9, 52.9, 64.0, 88.0, + 102.0, 154.0, 173.9, 195.0, 218.0, 240, 264, + 293, 317, 342, 367, 488, 520, 550, + 580, 640, 680, 760, 800, 850, 920, + 980, 1050, 1110, 1180, 1250, 1310, 1390, + 1460, 1520, 1590, 1660, 1720, 2033, 2094, + 2158, 2227, 2293, 2357, 2467, 2535, 2606, + 2674, 2944, 3010, 3082, 3149, 3433, 3510, + 3699, 3777, 6169, 6318, 6476, 6646, 6807, + 6964, 7283, 7450, 7630, 7800, 8260, 8410, + 8570, 8710, 9610, 9780, 10120, 10290, 21770, + 22160, 22600, 22990, 26310, 26830, 27490, 27903.1, 110842.0, 112843.7 }; #endif // #ifndef WARPX_IONIZATION_TABLE_H_ diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H index 8cb105aa0..a3faf95a0 100644 --- a/Source/Utils/NCIGodfreyTables.H +++ b/Source/Utils/NCIGodfreyTables.H @@ -12,212 +12,212 @@ const int tab_length = 101; // Table of coefficient for Ex, Ey abd Bz // We typically interpolate between two lines const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ - -2.47536,2.04288,-0.598163,0.0314711, - -2.47536,2.04288,-0.598163,0.0314711, - -2.47545,2.04309,-0.598307,0.0315029, - -2.4756,2.04342,-0.598549,0.0315558, - -2.47581,2.0439,-0.598886,0.0316298, - -2.47608,2.0445,-0.59932,0.031725, - -2.47641,2.04525,-0.59985,0.0318412, - -2.4768,2.04612,-0.600477,0.0319785, - -2.47725,2.04714,-0.6012,0.0321367, - -2.47776,2.04829,-0.602019,0.0323158, - -2.47833,2.04957,-0.602934,0.0325158, - -2.47896,2.05099,-0.603944,0.0327364, - -2.47965,2.05254,-0.605051,0.0329777, - -2.4804,2.05423,-0.606253,0.0332396, - -2.48121,2.05606,-0.60755,0.0335218, - -2.48208,2.05802,-0.608942,0.0338243, - -2.48301,2.06012,-0.610429,0.0341469, - -2.48401,2.06235,-0.61201,0.0344895, - -2.48506,2.06471,-0.613685,0.0348519, - -2.48618,2.06721,-0.615453,0.0352339, - -2.48735,2.06984,-0.617314,0.0356353, - -2.48859,2.07261,-0.619268,0.0360559, - -2.48988,2.0755,-0.621312,0.0364954, - -2.49123,2.07853,-0.623447,0.0369536, - -2.49265,2.08169,-0.625672,0.0374302, - -2.49412,2.08498,-0.627986,0.0379248, - -2.49565,2.0884,-0.630386,0.0384372, - -2.49724,2.09194,-0.632873,0.0389669, - -2.49888,2.09561,-0.635443,0.0395135, - -2.50058,2.09939,-0.638096,0.0400766, - -2.50234,2.1033,-0.640829,0.0406557, - -2.50415,2.10732,-0.64364,0.0412502, - -2.50601,2.11145,-0.646526,0.0418594, - -2.50791,2.1157,-0.649485,0.0424828, - -2.50987,2.12004,-0.652512,0.0431196, - -2.51187,2.12448,-0.655604,0.0437688, - -2.51392,2.12901,-0.658756,0.0444297, - -2.516,2.13363,-0.661964,0.0451011, - -2.51812,2.13832,-0.665221,0.0457818, - -2.52027,2.14308,-0.668521,0.0464705, - -2.52244,2.14789,-0.671856,0.0471658, - -2.52464,2.15274,-0.675218,0.0478658, - -2.52684,2.15762,-0.678596,0.0485687, - -2.52906,2.16251,-0.68198,0.0492723, - -2.53126,2.16738,-0.685355,0.049974, - -2.53345,2.17222,-0.688706,0.0506708, - -2.53561,2.177,-0.692015,0.0513594, - -2.53773,2.18168,-0.69526,0.0520359, - -2.53978,2.18623,-0.698416,0.0526955, - -2.54175,2.19059,-0.701452,0.053333, - -2.5436,2.19471,-0.704331,0.0539417, - -2.54531,2.19852,-0.70701,0.0545141, - -2.54683,2.20193,-0.709433,0.0550409, - -2.5481,2.20483,-0.711533,0.0555106, - -2.54906,2.20709,-0.713224,0.0559094, - -2.54963,2.20852,-0.714397,0.0562198, - -2.54968,2.20888,-0.714907,0.0564196, - -2.54905,2.20785,-0.714562,0.0564797, - -2.54751,2.20496,-0.713094,0.0563618, - -2.54472,2.19955,-0.710118,0.0560124, - -2.54014,2.19058,-0.705048,0.0553544, - -2.53286,2.1763,-0.69693,0.0542684, - -2.52115,2.15344,-0.684027,0.05255, - -2.50098,2.11466,-0.66255,0.0497817, - -2.45797,2.03459,-0.620099,0.0446889, - -2.28371,1.72254,-0.465905,0.0283268, - -2.4885,2.04899,-0.599292,0.0390466, - -2.1433,1.36735,-0.220924,-0.00215633, - -2.4943,2.07019,-0.610552,0.035166, - -2.84529,2.77303,-1.00018,0.0724884, - -2.72242,2.51888,-0.847226,0.0509964, - -2.65633,2.3744,-0.750392,0.0326366, - -2.59601,2.23412,-0.646421,0.00868027, - -2.51477,2.0369,-0.491066,-0.0306397, - -2.35935,1.65155,-0.178971,-0.112713, - -1.84315,0.361693,0.876104,-0.393844, - -2.65422,2.39262,-0.789663,0.0516265, - -3.46529,4.42354,-2.45543,0.497097, - -3.15747,3.65311,-1.824,0.328432, - -3.04694,3.37613,-1.59668,0.267631, - -2.99205,3.23814,-1.48302,0.237103, - -2.96075,3.15894,-1.41733,0.219317, - -2.94172,3.11028,-1.37649,0.20811, - -2.92994,3.07962,-1.35025,0.200755, - -2.92283,3.06054,-1.33338,0.195859, - -2.91894,3.04938,-1.3229,0.192637, - -2.91736,3.04394,-1.31702,0.190612, - -2.91753,3.04278,-1.31456,0.189477, - -2.91905,3.04494,-1.31475,0.189026, - -2.92165,3.04973,-1.31705,0.189117, - -2.92512,3.05667,-1.32105,0.189646, - -2.92933,3.06539,-1.32646,0.190538, - -2.93416,3.07562,-1.33308,0.191735, - -2.93952,3.08715,-1.34072,0.193194, - -2.94535,3.09982,-1.34925,0.194881, - -2.95159,3.11349,-1.35858,0.196769, - -2.9582,3.12805,-1.36861,0.198838, - -2.96514,3.14342,-1.37929,0.201068, - -2.97239,3.15953,-1.39055,0.203448, - -2.97991,3.17632,-1.40234,0.205964, + -2.47536,2.04288,-0.598163,0.0314711, + -2.47536,2.04288,-0.598163,0.0314711, + -2.47545,2.04309,-0.598307,0.0315029, + -2.4756,2.04342,-0.598549,0.0315558, + -2.47581,2.0439,-0.598886,0.0316298, + -2.47608,2.0445,-0.59932,0.031725, + -2.47641,2.04525,-0.59985,0.0318412, + -2.4768,2.04612,-0.600477,0.0319785, + -2.47725,2.04714,-0.6012,0.0321367, + -2.47776,2.04829,-0.602019,0.0323158, + -2.47833,2.04957,-0.602934,0.0325158, + -2.47896,2.05099,-0.603944,0.0327364, + -2.47965,2.05254,-0.605051,0.0329777, + -2.4804,2.05423,-0.606253,0.0332396, + -2.48121,2.05606,-0.60755,0.0335218, + -2.48208,2.05802,-0.608942,0.0338243, + -2.48301,2.06012,-0.610429,0.0341469, + -2.48401,2.06235,-0.61201,0.0344895, + -2.48506,2.06471,-0.613685,0.0348519, + -2.48618,2.06721,-0.615453,0.0352339, + -2.48735,2.06984,-0.617314,0.0356353, + -2.48859,2.07261,-0.619268,0.0360559, + -2.48988,2.0755,-0.621312,0.0364954, + -2.49123,2.07853,-0.623447,0.0369536, + -2.49265,2.08169,-0.625672,0.0374302, + -2.49412,2.08498,-0.627986,0.0379248, + -2.49565,2.0884,-0.630386,0.0384372, + -2.49724,2.09194,-0.632873,0.0389669, + -2.49888,2.09561,-0.635443,0.0395135, + -2.50058,2.09939,-0.638096,0.0400766, + -2.50234,2.1033,-0.640829,0.0406557, + -2.50415,2.10732,-0.64364,0.0412502, + -2.50601,2.11145,-0.646526,0.0418594, + -2.50791,2.1157,-0.649485,0.0424828, + -2.50987,2.12004,-0.652512,0.0431196, + -2.51187,2.12448,-0.655604,0.0437688, + -2.51392,2.12901,-0.658756,0.0444297, + -2.516,2.13363,-0.661964,0.0451011, + -2.51812,2.13832,-0.665221,0.0457818, + -2.52027,2.14308,-0.668521,0.0464705, + -2.52244,2.14789,-0.671856,0.0471658, + -2.52464,2.15274,-0.675218,0.0478658, + -2.52684,2.15762,-0.678596,0.0485687, + -2.52906,2.16251,-0.68198,0.0492723, + -2.53126,2.16738,-0.685355,0.049974, + -2.53345,2.17222,-0.688706,0.0506708, + -2.53561,2.177,-0.692015,0.0513594, + -2.53773,2.18168,-0.69526,0.0520359, + -2.53978,2.18623,-0.698416,0.0526955, + -2.54175,2.19059,-0.701452,0.053333, + -2.5436,2.19471,-0.704331,0.0539417, + -2.54531,2.19852,-0.70701,0.0545141, + -2.54683,2.20193,-0.709433,0.0550409, + -2.5481,2.20483,-0.711533,0.0555106, + -2.54906,2.20709,-0.713224,0.0559094, + -2.54963,2.20852,-0.714397,0.0562198, + -2.54968,2.20888,-0.714907,0.0564196, + -2.54905,2.20785,-0.714562,0.0564797, + -2.54751,2.20496,-0.713094,0.0563618, + -2.54472,2.19955,-0.710118,0.0560124, + -2.54014,2.19058,-0.705048,0.0553544, + -2.53286,2.1763,-0.69693,0.0542684, + -2.52115,2.15344,-0.684027,0.05255, + -2.50098,2.11466,-0.66255,0.0497817, + -2.45797,2.03459,-0.620099,0.0446889, + -2.28371,1.72254,-0.465905,0.0283268, + -2.4885,2.04899,-0.599292,0.0390466, + -2.1433,1.36735,-0.220924,-0.00215633, + -2.4943,2.07019,-0.610552,0.035166, + -2.84529,2.77303,-1.00018,0.0724884, + -2.72242,2.51888,-0.847226,0.0509964, + -2.65633,2.3744,-0.750392,0.0326366, + -2.59601,2.23412,-0.646421,0.00868027, + -2.51477,2.0369,-0.491066,-0.0306397, + -2.35935,1.65155,-0.178971,-0.112713, + -1.84315,0.361693,0.876104,-0.393844, + -2.65422,2.39262,-0.789663,0.0516265, + -3.46529,4.42354,-2.45543,0.497097, + -3.15747,3.65311,-1.824,0.328432, + -3.04694,3.37613,-1.59668,0.267631, + -2.99205,3.23814,-1.48302,0.237103, + -2.96075,3.15894,-1.41733,0.219317, + -2.94172,3.11028,-1.37649,0.20811, + -2.92994,3.07962,-1.35025,0.200755, + -2.92283,3.06054,-1.33338,0.195859, + -2.91894,3.04938,-1.3229,0.192637, + -2.91736,3.04394,-1.31702,0.190612, + -2.91753,3.04278,-1.31456,0.189477, + -2.91905,3.04494,-1.31475,0.189026, + -2.92165,3.04973,-1.31705,0.189117, + -2.92512,3.05667,-1.32105,0.189646, + -2.92933,3.06539,-1.32646,0.190538, + -2.93416,3.07562,-1.33308,0.191735, + -2.93952,3.08715,-1.34072,0.193194, + -2.94535,3.09982,-1.34925,0.194881, + -2.95159,3.11349,-1.35858,0.196769, + -2.9582,3.12805,-1.36861,0.198838, + -2.96514,3.14342,-1.37929,0.201068, + -2.97239,3.15953,-1.39055,0.203448, + -2.97991,3.17632,-1.40234,0.205964, -2.98769,3.19374,-1.41463,0.208607 }; // Table of coefficient for Bx, By and Ez // We typically interpolate between two lines const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ - -2.80862,2.80104,-1.14615,0.154077, - -2.80862,2.80104,-1.14615,0.154077, - -2.80851,2.80078,-1.14595,0.154027, - -2.80832,2.80034,-1.14561,0.153945, - -2.80807,2.79973,-1.14514,0.153829, - -2.80774,2.79894,-1.14454,0.15368, - -2.80733,2.79798,-1.1438,0.153498, - -2.80685,2.79685,-1.14292,0.153284, - -2.8063,2.79554,-1.14192,0.153036, - -2.80568,2.79405,-1.14077,0.152756, - -2.80498,2.79239,-1.1395,0.152443, - -2.80421,2.79056,-1.13809,0.152098, - -2.80337,2.78856,-1.13656,0.151721, - -2.80246,2.78638,-1.13488,0.151312, - -2.80147,2.78404,-1.13308,0.150871, - -2.80041,2.78152,-1.13115,0.150397, - -2.79927,2.77882,-1.12908,0.149893, - -2.79807,2.77596,-1.12689,0.149356, - -2.79679,2.77292,-1.12456,0.148789, - -2.79543,2.76972,-1.12211,0.14819, - -2.79401,2.76634,-1.11953,0.14756, - -2.79251,2.76279,-1.11681,0.1469, - -2.79094,2.75907,-1.11397,0.146208, - -2.78929,2.75517,-1.111,0.145486, - -2.78757,2.7511,-1.10789,0.144733, - -2.78578,2.74686,-1.10466,0.14395, - -2.78391,2.74245,-1.1013,0.143137, - -2.78196,2.73786,-1.09781,0.142293, - -2.77994,2.73309,-1.09419,0.141419, - -2.77784,2.72814,-1.09043,0.140514, - -2.77566,2.72301,-1.08654,0.139578, - -2.7734,2.7177,-1.08252,0.138612, - -2.77106,2.7122,-1.07836,0.137614, - -2.76864,2.70651,-1.07406,0.136586, - -2.76613,2.70062,-1.06962,0.135525, - -2.76353,2.69453,-1.06503,0.134432, - -2.76084,2.68824,-1.0603,0.133307, - -2.75806,2.68173,-1.05541,0.132148, - -2.75518,2.675,-1.05037,0.130954, - -2.75219,2.66804,-1.04516,0.129725, - -2.7491,2.66084,-1.03978,0.12846, - -2.7459,2.65339,-1.03423,0.127156, - -2.74257,2.64566,-1.02848,0.125813, - -2.73912,2.63765,-1.02254,0.124428, - -2.73552,2.62934,-1.01638,0.122999, - -2.73178,2.62069,-1.01,0.121523, - -2.72787,2.61169,-1.00337,0.119996, - -2.72379,2.6023,-0.996479,0.118417, - -2.71951,2.59248,-0.989294,0.116778, - -2.71501,2.58218,-0.981786,0.115076, - -2.71026,2.57135,-0.97392,0.113303, - -2.70524,2.55991,-0.965651,0.111453, - -2.69989,2.54778,-0.956922,0.109514, - -2.69416,2.53484,-0.947666,0.107476, - -2.68799,2.52096,-0.937795,0.105324, - -2.68129,2.50596,-0.927197,0.103039, - -2.67394,2.48959,-0.915724,0.100597, - -2.66578,2.47153,-0.903179,0.097968, - -2.65657,2.4513,-0.889283,0.0951084, - -2.64598,2.42824,-0.873638,0.0919592, - -2.63347,2.40127,-0.855632,0.0884325, - -2.61813,2.36864,-0.834261,0.0843898, - -2.59821,2.32701,-0.807691,0.0795876, - -2.56971,2.26887,-0.77188,0.0735132, - -2.51823,2.16823,-0.713448,0.0645399, - -2.33537,1.8294,-0.533852,0.0409941, - -2.53143,2.14818,-0.670502,0.053982, - -2.17737,1.43641,-0.259095,0.00101255, - -2.51929,2.12931,-0.654743,0.0452381, - -2.86122,2.82221,-1.05039,0.0894636, - -2.72908,2.54506,-0.87834,0.0626188, - -2.6536,2.37954,-0.7665,0.0409117, - -2.58374,2.21923,-0.649738,0.0146791, - -2.49284,2.00346,-0.48457,-0.0255348, - -2.32762,1.60337,-0.1698,-0.105287, - -1.80149,0.316787,0.855414,-0.369652, - -2.60242,2.28418,-0.721378,0.040091, - -3.40335,4.25157,-2.29817,0.449834, - -3.0852,3.47341,-1.67791,0.28982, - -2.9642,3.17856,-1.44399,0.229852, - -2.89872,3.01966,-1.31861,0.197945, - -2.85668,2.91811,-1.23894,0.17783, - -2.82679,2.84621,-1.18287,0.163785, - -2.80401,2.79167,-1.14058,0.153278, - -2.78577,2.74819,-1.10706,0.145015, - -2.77061,2.7122,-1.07946,0.138267, - -2.75764,2.68152,-1.05606,0.132589, - -2.74627,2.65475,-1.03575,0.127695, - -2.73612,2.63093,-1.01777,0.123395, - -2.72692,2.6094,-1.00159,0.119553, - -2.71846,2.58968,-0.986841,0.116074, - -2.71061,2.57142,-0.973239,0.112887, - -2.70323,2.55434,-0.960573,0.109937, - -2.69626,2.53824,-0.948678,0.107185, - -2.68962,2.52294,-0.937429,0.104598, - -2.68327,2.50833,-0.926722,0.102151, - -2.67714,2.4943,-0.916477,0.0998223, - -2.67122,2.48076,-0.906627,0.0975966, - -2.66546,2.46764,-0.897118,0.0954599, - -2.65985,2.45489,-0.887903,0.0934011, + -2.80862,2.80104,-1.14615,0.154077, + -2.80862,2.80104,-1.14615,0.154077, + -2.80851,2.80078,-1.14595,0.154027, + -2.80832,2.80034,-1.14561,0.153945, + -2.80807,2.79973,-1.14514,0.153829, + -2.80774,2.79894,-1.14454,0.15368, + -2.80733,2.79798,-1.1438,0.153498, + -2.80685,2.79685,-1.14292,0.153284, + -2.8063,2.79554,-1.14192,0.153036, + -2.80568,2.79405,-1.14077,0.152756, + -2.80498,2.79239,-1.1395,0.152443, + -2.80421,2.79056,-1.13809,0.152098, + -2.80337,2.78856,-1.13656,0.151721, + -2.80246,2.78638,-1.13488,0.151312, + -2.80147,2.78404,-1.13308,0.150871, + -2.80041,2.78152,-1.13115,0.150397, + -2.79927,2.77882,-1.12908,0.149893, + -2.79807,2.77596,-1.12689,0.149356, + -2.79679,2.77292,-1.12456,0.148789, + -2.79543,2.76972,-1.12211,0.14819, + -2.79401,2.76634,-1.11953,0.14756, + -2.79251,2.76279,-1.11681,0.1469, + -2.79094,2.75907,-1.11397,0.146208, + -2.78929,2.75517,-1.111,0.145486, + -2.78757,2.7511,-1.10789,0.144733, + -2.78578,2.74686,-1.10466,0.14395, + -2.78391,2.74245,-1.1013,0.143137, + -2.78196,2.73786,-1.09781,0.142293, + -2.77994,2.73309,-1.09419,0.141419, + -2.77784,2.72814,-1.09043,0.140514, + -2.77566,2.72301,-1.08654,0.139578, + -2.7734,2.7177,-1.08252,0.138612, + -2.77106,2.7122,-1.07836,0.137614, + -2.76864,2.70651,-1.07406,0.136586, + -2.76613,2.70062,-1.06962,0.135525, + -2.76353,2.69453,-1.06503,0.134432, + -2.76084,2.68824,-1.0603,0.133307, + -2.75806,2.68173,-1.05541,0.132148, + -2.75518,2.675,-1.05037,0.130954, + -2.75219,2.66804,-1.04516,0.129725, + -2.7491,2.66084,-1.03978,0.12846, + -2.7459,2.65339,-1.03423,0.127156, + -2.74257,2.64566,-1.02848,0.125813, + -2.73912,2.63765,-1.02254,0.124428, + -2.73552,2.62934,-1.01638,0.122999, + -2.73178,2.62069,-1.01,0.121523, + -2.72787,2.61169,-1.00337,0.119996, + -2.72379,2.6023,-0.996479,0.118417, + -2.71951,2.59248,-0.989294,0.116778, + -2.71501,2.58218,-0.981786,0.115076, + -2.71026,2.57135,-0.97392,0.113303, + -2.70524,2.55991,-0.965651,0.111453, + -2.69989,2.54778,-0.956922,0.109514, + -2.69416,2.53484,-0.947666,0.107476, + -2.68799,2.52096,-0.937795,0.105324, + -2.68129,2.50596,-0.927197,0.103039, + -2.67394,2.48959,-0.915724,0.100597, + -2.66578,2.47153,-0.903179,0.097968, + -2.65657,2.4513,-0.889283,0.0951084, + -2.64598,2.42824,-0.873638,0.0919592, + -2.63347,2.40127,-0.855632,0.0884325, + -2.61813,2.36864,-0.834261,0.0843898, + -2.59821,2.32701,-0.807691,0.0795876, + -2.56971,2.26887,-0.77188,0.0735132, + -2.51823,2.16823,-0.713448,0.0645399, + -2.33537,1.8294,-0.533852,0.0409941, + -2.53143,2.14818,-0.670502,0.053982, + -2.17737,1.43641,-0.259095,0.00101255, + -2.51929,2.12931,-0.654743,0.0452381, + -2.86122,2.82221,-1.05039,0.0894636, + -2.72908,2.54506,-0.87834,0.0626188, + -2.6536,2.37954,-0.7665,0.0409117, + -2.58374,2.21923,-0.649738,0.0146791, + -2.49284,2.00346,-0.48457,-0.0255348, + -2.32762,1.60337,-0.1698,-0.105287, + -1.80149,0.316787,0.855414,-0.369652, + -2.60242,2.28418,-0.721378,0.040091, + -3.40335,4.25157,-2.29817,0.449834, + -3.0852,3.47341,-1.67791,0.28982, + -2.9642,3.17856,-1.44399,0.229852, + -2.89872,3.01966,-1.31861,0.197945, + -2.85668,2.91811,-1.23894,0.17783, + -2.82679,2.84621,-1.18287,0.163785, + -2.80401,2.79167,-1.14058,0.153278, + -2.78577,2.74819,-1.10706,0.145015, + -2.77061,2.7122,-1.07946,0.138267, + -2.75764,2.68152,-1.05606,0.132589, + -2.74627,2.65475,-1.03575,0.127695, + -2.73612,2.63093,-1.01777,0.123395, + -2.72692,2.6094,-1.00159,0.119553, + -2.71846,2.58968,-0.986841,0.116074, + -2.71061,2.57142,-0.973239,0.112887, + -2.70323,2.55434,-0.960573,0.109937, + -2.69626,2.53824,-0.948678,0.107185, + -2.68962,2.52294,-0.937429,0.104598, + -2.68327,2.50833,-0.926722,0.102151, + -2.67714,2.4943,-0.916477,0.0998223, + -2.67122,2.48076,-0.906627,0.0975966, + -2.66546,2.46764,-0.897118,0.0954599, + -2.65985,2.45489,-0.887903,0.0934011, -2.65437,2.44244,-0.878945,0.0914107 }; diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index a5ab7b18c..f84701a02 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -76,7 +76,7 @@ WarpX::MoveWindow (bool move_j) new_slice_lo[i] = current_slice_lo[i]; new_slice_hi[i] = current_slice_hi[i]; } - int num_shift_base_slice = static_cast<int> ((moving_window_x - + int num_shift_base_slice = static_cast<int> ((moving_window_x - current_slice_lo[dir]) / cdx[dir]); new_slice_lo[dir] = current_slice_lo[dir] + num_shift_base_slice*cdx[dir]; new_slice_hi[dir] = current_slice_hi[dir] + num_shift_base_slice*cdx[dir]; diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index 39fded8e8..fd6e72dc6 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -7,5 +7,5 @@ void ReadBoostedFrameParameters(amrex::Real& gamma_boost, amrex::Real& beta_boos void ConvertLabParamsToBoost(); -void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, +void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax); diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 19e898208..4b11eb69d 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -68,7 +68,7 @@ void ConvertLabParamsToBoost() BL_ASSERT(slice_lo.size() == AMREX_SPACEDIM); pp_slice.queryarr("dom_hi",slice_hi,0,AMREX_SPACEDIM); BL_ASSERT(slice_hi.size() == AMREX_SPACEDIM); - + pp_amr.query("max_level", max_level); if (max_level > 0){ @@ -113,7 +113,7 @@ void ConvertLabParamsToBoost() } -/* \brief Function that sets the value of MultiFab MF to zero for z between +/* \brief Function that sets the value of MultiFab MF to zero for z between * zmin and zmax. */ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax){ @@ -143,7 +143,7 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) const Real z_gridpoint = zmin_box+(k-lo_ind)*dz; #else const Real z_gridpoint = zmin_box+(j-lo_ind)*dz; -#endif +#endif if ( (z_gridpoint >= zmin) && (z_gridpoint < zmax) ) { arr(i,j,k) = 0.; }; diff --git a/Source/Utils/utils_ES.F90 b/Source/Utils/utils_ES.F90 index 887a5ef15..ce143bb94 100644 --- a/Source/Utils/utils_ES.F90 +++ b/Source/Utils/utils_ES.F90 @@ -16,9 +16,9 @@ contains integer, intent(in) :: lrat(3) real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2),clo(3):chi(3)) real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2),flo(3):fhi(3)) - + integer :: i, j, k, ii, jj, kk - + do k = lo(3), hi(3) kk = k * lrat(3) do j = lo(2), hi(2) @@ -27,7 +27,7 @@ contains ii = i * lrat(1) crse(i,j,k) = fine(ii,jj,kk) + & ! These six fine nodes are shared by two coarse nodes... - 0.5d0 * (fine(ii-1,jj,kk) + fine(ii+1,jj,kk) + & + 0.5d0 * (fine(ii-1,jj,kk) + fine(ii+1,jj,kk) + & fine(ii,jj-1,kk) + fine(ii,jj+1,kk) + & fine(ii,jj,kk-1) + fine(ii,jj,kk+1)) + & ! ... these twelve are shared by four... @@ -59,16 +59,16 @@ contains integer, intent(in) :: lrat(2) real(amrex_real), intent(inout) :: crse(clo(1):chi(1),clo(2):chi(2)) real(amrex_real), intent(in) :: fine(flo(1):fhi(1),flo(2):fhi(2)) - + integer :: i, j, ii, jj - + do j = lo(2), hi(2) jj = j * lrat(2) do i = lo(1), hi(1) ii = i * lrat(1) crse(i,j) = fine(ii,jj) + & ! These four fine nodes are shared by two coarse nodes... - 0.5d0 * (fine(ii-1,jj) + fine(ii+1,jj) + & + 0.5d0 * (fine(ii-1,jj) + fine(ii+1,jj) + & fine(ii,jj-1) + fine(ii,jj+1)) + & ! ... and these four are shared by four... 0.25d0 * (fine(ii-1,jj-1) + fine(ii-1,jj+1) + & @@ -87,12 +87,12 @@ contains double precision, intent(inout) :: input_data(lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) double precision, intent(inout) :: bndry_data(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1,lo(3)-1:hi(3)+1) integer(c_int), intent(in ) :: mask (lo(1):hi(1),lo(2):hi(2),lo(3):hi(3)) - + integer :: i, j, k - + do k = lo(3), hi(3) do j = lo(2), hi(2) - do i = lo(1), hi(1) + do i = lo(1), hi(1) if (mask(i,j,k) .eq. 1) then bndry_data(i,j,k) = input_data(i,j,k) input_data(i,j,k) = 0.d0 @@ -110,11 +110,11 @@ contains double precision, intent(inout) :: input_data(lo(1):hi(1),lo(2):hi(2)) double precision, intent(inout) :: bndry_data(lo(1)-1:hi(1)+1,lo(2)-1:hi(2)+1) integer(c_int), intent(in ) :: mask (lo(1):hi(1),lo(2):hi(2)) - + integer :: i, j - + do j = lo(2), hi(2) - do i = lo(1), hi(1) + do i = lo(1), hi(1) if (mask(i,j) .eq. 1) then bndry_data(i,j) = input_data(i,j) input_data(i,j) = 0.d0 @@ -145,7 +145,7 @@ contains end do end do end do - + if (total .gt. 0) then mask(i,j,k) = 1 else @@ -169,14 +169,14 @@ contains do j = lo(2), hi(2) do i = lo(1), hi(1) - + total = 0 do ii = i-ncells, i+ncells do jj = j-ncells, j+ncells total = total + tmp_mask(ii, jj) end do end do - + if (total .gt. 0) then mask(i,j) = 1 else diff --git a/Source/Utils/write_atomic_data_cpp.py b/Source/Utils/write_atomic_data_cpp.py index b085d50eb..1c85fa32c 100644 --- a/Source/Utils/write_atomic_data_cpp.py +++ b/Source/Utils/write_atomic_data_cpp.py @@ -60,7 +60,7 @@ for element in ion_names: '\n\s+(\d+)\s+\|\s+%s\s+\w+\s+\|\s+\+*(\d+)\s+\|\s+\(*\[*(\d+\.*\d*)' \ %element list_of_tuples = re.findall( regex_command, text_data ) - for count, energy in enumerate([x[2] for x in list_of_tuples]): + for count, energy in enumerate([x[2] for x in list_of_tuples]): if count%7==0: cpp_string += '\n ' cpp_string += energy + ', ' cpp_string = cpp_string[:-2] diff --git a/Source/main.cpp b/Source/main.cpp index d31268b12..6aca49cd3 100644 --- a/Source/main.cpp +++ b/Source/main.cpp @@ -25,18 +25,18 @@ int main(int argc, char* argv[]) ConvertLabParamsToBoost(); BL_PROFILE_VAR("main()", pmain); - + const Real strt_total = amrex::second(); { WarpX warpx; - + warpx.InitData(); warpx.Evolve(); - + Real end_total = amrex::second() - strt_total; - + ParallelDescriptor::ReduceRealMax(end_total ,ParallelDescriptor::IOProcessorNumber()); if (warpx.Verbose()) { amrex::Print() << "Total Time : " << end_total << '\n'; |