diff options
Diffstat (limited to 'Source')
60 files changed, 3097 insertions, 985 deletions
diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index c3449cecd..f780f335c 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -412,7 +412,7 @@ PML::MakeBoxArray (const amrex::Geometry& geom, const amrex::BoxArray& grid_ba, { Box domain = geom.Domain(); for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { - if ( ! Geometry::isPeriodic(idim) ) { + if ( ! geom.isPeriodic(idim) ) { domain.grow(idim, ncell); } } diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index e35d307a6..3a09259b0 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -2,11 +2,13 @@ #define WARPX_BoostedFrameDiagnostic_H_ #include <vector> +#include <map> #include <AMReX_VisMF.H> #include <AMReX_PlotFileUtil.H> #include <AMReX_ParallelDescriptor.H> #include <AMReX_Geometry.H> +#include <AMReX_CudaContainers.H> #include "MultiParticleContainer.H" #include "WarpXConst.H" @@ -32,17 +34,25 @@ class BoostedFrameDiagnostic { std::string file_name; amrex::Real t_lab; - amrex::Real zmin_lab; - amrex::Real zmax_lab; + amrex::RealBox prob_domain_lab_; + amrex::IntVect prob_ncells_lab_; + amrex::Real current_z_lab; amrex::Real current_z_boost; + + int ncomp_to_dump_; + std::vector<std::string> mesh_field_names_; + int file_num; int initial_i; const BoostedFrameDiagnostic& my_bfd; LabSnapShot(amrex::Real t_lab_in, amrex::Real t_boost, - amrex::Real zmin_lab_in, - amrex::Real zmax_lab_in, int file_num_in, + const amrex::RealBox prob_domain_lab, + const amrex::IntVect prob_ncells_lab, + const int ncomp_to_dump, + const std::vector<std::string> mesh_field_names, + int file_num_in, const BoostedFrameDiagnostic& bfd); /// @@ -69,15 +79,21 @@ class BoostedFrameDiagnostic { amrex::Real dt_snapshots_lab_; amrex::Real dt_boost_; int N_snapshots_; - unsigned Nx_lab_; - unsigned Ny_lab_; - unsigned Nz_lab_; int boost_direction_; + // For back-transformed diagnostics of grid fields, data_buffer_[i] + // stores a buffer of the fields in the lab frame (in a MultiFab, i.e. + // with all box data etc.). When the buffer if full, dump to file. amrex::Vector<std::unique_ptr<amrex::MultiFab> > data_buffer_; + // particles_buffer_ is currently blind to refinement level. + // particles_buffer_[i][j] is a WarpXParticleContainer::DiagnosticParticleData where + // - i is the back-transformed snapshot number + // - j is the species number amrex::Vector<amrex::Vector<WarpXParticleContainer::DiagnosticParticleData> > particles_buffer_; int num_buffer_ = 256; int max_box_size_ = 256; + // buff_counter_[i] is the number of z slices in data_buffer_[i] + // for snapshot number i. amrex::Vector<int> buff_counter_; amrex::Vector<LabSnapShot> snapshots_; @@ -105,6 +121,32 @@ public: 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}, + {"rho", 9} }; + + // 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", + "jx", "jy", "jz", "rho"}; + int ncomp_to_dump = 10; + }; #endif diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 13972075d..709b7cb48 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -17,8 +17,6 @@ using namespace amrex; namespace { const std::vector<std::string> particle_field_names = {"w", "x", "y", "z", "ux", "uy", "uz"}; - const std::vector<std::string> mesh_field_names = - {"Ex", "Ey", "Ez", "Bx", "By", "Bz", "jx", "jy", "jz", "rho"}; /* Creates the HDF5 file in truncate mode and closes it. @@ -446,6 +444,83 @@ namespace } #endif +namespace +{ + void + CopySlice(MultiFab& tmp, MultiFab& buf, int k_lab, + const Gpu::ManagedDeviceVector<int>& map_actual_fields_to_dump) + { + const int ncomp_to_dump = map_actual_fields_to_dump.size(); + // Copy data from MultiFab tmp to MultiFab data_buffer[i]. +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(tmp, TilingIfNotGPU()); mfi.isValid(); ++mfi) { + Array4< Real> tmp_arr = tmp[mfi].array(); + Array4< Real> buf_arr = buf[mfi].array(); + // For 3D runs, tmp is a 2D (x,y) multifab, that contains only + // slice to write to file. + const Box& bx = mfi.tilebox(); + + const auto field_map_ptr = map_actual_fields_to_dump.dataPtr(); + ParallelFor(bx, ncomp_to_dump, + [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) + { + const int icomp = field_map_ptr[n]; +#if (AMREX_SPACEDIM == 3) + buf_arr(i,j,k_lab,n) += tmp_arr(i,j,k,icomp); +#else + buf_arr(i,k_lab,k,n) += tmp_arr(i,j,k,icomp); +#endif + } + ); + } + } + +void +LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) +{ + // Loop over tiles/boxes and in-place convert each slice from boosted + // frame to lab frame. +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + 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 + // 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 + // changed by the tranform. + Real e_lab, b_lab, j_lab, r_lab; + e_lab = gamma_boost * (arr(i, j, k, 0) + beta_boost*clight*arr(i, j, k, 4)); + b_lab = gamma_boost * (arr(i, j, k, 4) + beta_boost*arr(i, j, k, 0)/clight); + + arr(i, j, k, 0) = e_lab; + arr(i, j, k, 4) = b_lab; + + e_lab = gamma_boost * (arr(i, j, k, 1) - beta_boost*clight*arr(i, j, k, 3)); + b_lab = gamma_boost * (arr(i, j, k, 3) - beta_boost*arr(i, j, k, 1)/clight); + + arr(i, j, k, 1) = e_lab; + arr(i, j, k, 3) = b_lab; + + // Transform the charge and current density. Only the z component of j is affected. + j_lab = gamma_boost*(arr(i, j, k, 8) + beta_boost*clight*arr(i, j, k, 9)); + r_lab = gamma_boost*(arr(i, j, k, 9) + beta_boost*arr(i, j, k, 8)/clight); + + arr(i, j, k, 8) = j_lab; + arr(i, j, k, 9) = r_lab; + } + ); + } +} +} + BoostedFrameDiagnostic:: BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, Real dt_snapshots_lab, int N_snapshots, @@ -469,23 +544,53 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, dz_lab_ = PhysConst::c * dt_boost_ * inv_beta_boost_ * inv_gamma_boost_; inv_dz_lab_ = 1.0 / dz_lab_; - Nx_lab_ = geom.Domain().length(0); + int Nz_lab = static_cast<unsigned>((zmax_lab - zmin_lab) * inv_dz_lab_); + int Nx_lab = geom.Domain().length(0); #if (AMREX_SPACEDIM == 3) - Ny_lab_ = geom.Domain().length(1); + int Ny_lab = geom.Domain().length(1); + IntVect prob_ncells_lab = {Nx_lab, Ny_lab, Nz_lab}; #else - Ny_lab_ = 1; + // Ny_lab = 1; + IntVect prob_ncells_lab = {Nx_lab, Nz_lab}; #endif - Nz_lab_ = static_cast<unsigned>((zmax_lab - zmin_lab) * inv_dz_lab_); - + writeMetaData(); if (WarpX::do_boosted_frame_fields) data_buffer_.resize(N_snapshots); if (WarpX::do_boosted_frame_particles) particles_buffer_.resize(N_snapshots); + + // Query fields to dump + std::vector<std::string> user_fields_to_dump; + ParmParse pp("warpx"); + bool do_user_fields; + do_user_fields = pp.queryarr("boosted_frame_diag_fields", + user_fields_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){ + ncomp_to_dump = user_fields_to_dump.size(); + map_actual_fields_to_dump.resize(ncomp_to_dump); + mesh_field_names.resize(ncomp_to_dump); + for (int i=0; i<ncomp_to_dump; i++){ + std::string fieldstr = user_fields_to_dump[i]; + mesh_field_names[i] = fieldstr; + map_actual_fields_to_dump[i] = possible_fields_to_dump[fieldstr]; + } + } + for (int i = 0; i < N_snapshots; ++i) { Real t_lab = i * dt_snapshots_lab_; - LabSnapShot snapshot(t_lab, t_boost, - zmin_lab + v_window_lab * t_lab, - zmax_lab + v_window_lab * t_lab, i, *this); + // Get simulation domain physical coordinates (in boosted frame). + RealBox prob_domain_lab = geom.ProbDomain(); + // Replace z bounds by lab-frame coordinates + // x and y bounds are the same for lab frame and boosted frame + prob_domain_lab.setLo(AMREX_SPACEDIM-1, zmin_lab + v_window_lab * t_lab); + prob_domain_lab.setHi(AMREX_SPACEDIM-1, zmax_lab + v_window_lab * t_lab); + // Construct LabSnapShot + LabSnapShot snapshot(t_lab, t_boost, prob_domain_lab, + prob_ncells_lab, ncomp_to_dump, + mesh_field_names, i, *this); snapshots_.push_back(snapshot); buff_counter_.push_back(0); if (WarpX::do_boosted_frame_fields) data_buffer_[i].reset( nullptr ); @@ -504,9 +609,11 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) 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) { - int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; + 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) { @@ -539,14 +646,20 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) } if (WarpX::do_boosted_frame_particles) { - for (int j = 0; j < mypc.nSpecies(); ++j) { + // Loop over species to be dumped to BFD + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + // Get species name + std::string species_name = + species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; #ifdef WARPX_USE_HDF5 + // Dump species data writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; + // Dump species data writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -575,73 +688,79 @@ writeLabFrameData(const MultiFab* cell_centered_data, const Real zhi_boost = domain_z_boost.hi(boost_direction_); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); - + + // 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, inv_gamma_boost_, inv_beta_boost_); + + Real zmin_lab = snapshots_[i].prob_domain_lab_.lo(AMREX_SPACEDIM-1); + Real zmax_lab = snapshots_[i].prob_domain_lab_.hi(AMREX_SPACEDIM-1); + // 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 < snapshots_[i].zmin_lab) or - (snapshots_[i].current_z_lab > snapshots_[i].zmax_lab) ) continue; + (snapshots_[i].current_z_lab < zmin_lab) or + (snapshots_[i].current_z_lab > zmax_lab) ) continue; - int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_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_; + // If buffer of snapshot i is empty... if (buff_counter_[i] == 0) { + // ... reset fields buffer data_buffer_[i] if (WarpX::do_boosted_frame_fields) { - const int ncomp = cell_centered_data->nComp(); Box buff_box = geom.Domain(); buff_box.setSmall(boost_direction_, i_lab - num_buffer_ + 1); buff_box.setBig(boost_direction_, i_lab); BoxArray buff_ba(buff_box); buff_ba.maxSize(max_box_size_); DistributionMapping buff_dm(buff_ba); - data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); + data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp_to_dump, 0) ); } - if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nSpecies()); + // ... reset particle buffer particles_buffer_[i] + if (WarpX::do_boosted_frame_particles) + particles_buffer_[i].resize(mypc.nSpeciesBoostedFrameDiags()); } if (WarpX::do_boosted_frame_fields) { const int ncomp = cell_centered_data->nComp(); const int start_comp = 0; const bool interpolate = true; + // Get slice in the boosted frame std::unique_ptr<MultiFab> slice = amrex::get_slice_data(boost_direction_, snapshots_[i].current_z_boost, *cell_centered_data, geom, start_comp, ncomp, interpolate); // transform it to the lab frame - for (MFIter mfi(*slice); mfi.isValid(); ++mfi) { - const Box& tile_box = mfi.tilebox(); - WRPX_LORENTZ_TRANSFORM_Z(BL_TO_FORTRAN_ANYD((*slice)[mfi]), - BL_TO_FORTRAN_BOX(tile_box), - &gamma_boost_, &beta_boost_); - } - + LorentzTransformZ(*slice, gamma_boost_, beta_boost_, ncomp); + // Create a 2D box for the slice in the boosted frame Real dx = geom.CellSize(boost_direction_); int i_boost = (snapshots_[i].current_z_boost - geom.ProbLo(boost_direction_))/dx; Box slice_box = geom.Domain(); slice_box.setSmall(boost_direction_, i_boost); slice_box.setBig(boost_direction_, i_boost); + // Make it a BoxArray slice_ba BoxArray slice_ba(slice_box); slice_ba.maxSize(max_box_size_); + // Create MultiFab tmp on slice_ba with data from slice MultiFab tmp(slice_ba, data_buffer_[i]->DistributionMap(), ncomp, 0); tmp.copy(*slice, 0, 0, ncomp); #ifdef _OPENMP #pragma omp parallel #endif - for (MFIter mfi(tmp, true); mfi.isValid(); ++mfi) { - const Box& tile_box = mfi.tilebox(); - WRPX_COPY_SLICE(BL_TO_FORTRAN_BOX(tile_box), - BL_TO_FORTRAN_ANYD(tmp[mfi]), - BL_TO_FORTRAN_ANYD((*data_buffer_[i])[mfi]), - &ncomp, &i_boost, &i_lab); - } + // Copy data from MultiFab tmp to MultiDab data_buffer[i] + CopySlice(tmp, *data_buffer_[i], i_lab, map_actual_fields_to_dump); } - + if (WarpX::do_boosted_frame_particles) { mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, old_z_boost, snapshots_[i].current_z_boost, @@ -651,6 +770,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, ++buff_counter_[i]; + // If buffer full, write to disk. if (buff_counter_[i] == num_buffer_) { if (WarpX::do_boosted_frame_fields) { @@ -666,14 +786,21 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) { - for (int j = 0; j < mypc.nSpecies(); ++j) { + // Loop over species to be dumped to BFD + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + // Get species name + const std::string species_name = species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; #ifdef WARPX_USE_HDF5 + // Write data to disk (HDF5) writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; + + // Write data to disk (custom) writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -791,15 +918,14 @@ writeMetaData () BL_PROFILE("BoostedFrameDiagnostic::writeMetaData"); if (ParallelDescriptor::IOProcessor()) { - std::string DiagnosticDirectory = "lab_frame_data"; - if (!UtilCreateDirectory(DiagnosticDirectory, 0755)) - CreateDirectoryFailed(DiagnosticDirectory); + if (!UtilCreateDirectory(WarpX::lab_data_directory, 0755)) + CreateDirectoryFailed(WarpX::lab_data_directory); VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size); std::ofstream HeaderFile; HeaderFile.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size()); - std::string HeaderFileName(DiagnosticDirectory + "/Header"); + std::string HeaderFileName(WarpX::lab_data_directory + "/Header"); HeaderFile.open(HeaderFileName.c_str(), std::ofstream::out | std::ofstream::trunc | std::ofstream::binary); @@ -812,25 +938,30 @@ writeMetaData () HeaderFile << dt_snapshots_lab_ << "\n"; HeaderFile << gamma_boost_ << "\n"; HeaderFile << beta_boost_ << "\n"; - HeaderFile << dz_lab_ << "\n"; - HeaderFile << Nz_lab_ << "\n"; } } BoostedFrameDiagnostic::LabSnapShot:: -LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, - Real zmax_lab_in, int file_num_in, const BoostedFrameDiagnostic& bfd) +LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, + IntVect prob_ncells_lab, + int ncomp_to_dump, + std::vector<std::string> mesh_field_names, + int file_num_in, + const BoostedFrameDiagnostic& bfd) : t_lab(t_lab_in), - zmin_lab(zmin_lab_in), - zmax_lab(zmax_lab_in), + prob_domain_lab_(prob_domain_lab), + prob_ncells_lab_(prob_ncells_lab), + ncomp_to_dump_(ncomp_to_dump), + mesh_field_names_(mesh_field_names), file_num(file_num_in), my_bfd(bfd) { + Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1); current_z_lab = 0.0; current_z_boost = 0.0; updateCurrentZPositions(t_boost, my_bfd.inv_gamma_boost_, my_bfd.inv_beta_boost_); initial_i = (current_z_lab - zmin_lab) / my_bfd.dz_lab_; - file_name = Concatenate("lab_frame_data/snapshot", file_num, 5); + file_name = Concatenate(WarpX::lab_data_directory + "/snapshot", file_num, 5); #ifdef WARPX_USE_HDF5 if (ParallelDescriptor::IOProcessor()) @@ -844,27 +975,34 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, { if (WarpX::do_boosted_frame_fields) { - for (int comp = 0; comp < static_cast<int>(mesh_field_names.size()); ++comp) { - output_create_field(file_name, mesh_field_names[comp], - my_bfd.Nx_lab_, - my_bfd.Ny_lab_, - my_bfd.Nz_lab_+1); + for (int comp = 0; comp < ncomp_to_dump; ++comp) { + output_create_field(file_name, mesh_field_names_[comp], + prob_ncells_lab_[0], +#if ( AMREX_SPACEDIM == 3 ) + prob_ncells_lab_[1], +#else + 1, +#endif + prob_ncells_lab_[AMREX_SPACEDIM-1]+1); } } } ParallelDescriptor::Barrier(); - if (WarpX::do_boosted_frame_particles) - { + if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); - for (int j = 0; j < mypc.nSpecies(); ++j) + // Loop over species to be dumped to BFD + for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { - output_create_species_group(file_name, species_names[j]); + // Loop over species to be dumped to BFD + 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) { - std::string field_path = species_names[j] + "/" + particle_field_names[k]; + std::string field_path = species_name + "/" + particle_field_names[k]; output_create_particle_field(file_name, field_path); } } @@ -884,11 +1022,14 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); - int nspecies = mypc.nSpecies(); const std::string particles_prefix = "particle"; - for(int i = 0; i < nspecies; ++i) { - const std::string fullpath = file_name + "/" + species_names[i]; + // Loop over species to be dumped to BFD + for(int i = 0; i < mypc.nSpeciesBoostedFrameDiags(); ++i) { + // Get species name + std::string species_name = + species_names[mypc.mapSpeciesBoostedFrameDiags(i)]; + const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } @@ -924,8 +1065,31 @@ writeSnapShotHeader() { HeaderFile.precision(17); HeaderFile << t_lab << "\n"; - HeaderFile << zmin_lab << "\n"; - HeaderFile << zmax_lab << "\n"; + // Write domain number of cells + HeaderFile << prob_ncells_lab_[0] << ' ' +#if ( AMREX_SPACEDIM==3 ) + << prob_ncells_lab_[1] << ' ' +#endif + << prob_ncells_lab_[AMREX_SPACEDIM-1] <<'\n'; + // Write domain physical boundaries + // domain lower bound + HeaderFile << prob_domain_lab_.lo(0) << ' ' +#if ( AMREX_SPACEDIM==3 ) + << prob_domain_lab_.lo(1) << ' ' +#endif + << prob_domain_lab_.lo(AMREX_SPACEDIM-1) <<'\n'; + // domain higher bound + HeaderFile << prob_domain_lab_.hi(0) << ' ' +#if ( AMREX_SPACEDIM==3 ) + << prob_domain_lab_.hi(1) << ' ' +#endif + << prob_domain_lab_.hi(AMREX_SPACEDIM-1) <<'\n'; + // List of fields dumped to file + for (int i=0; i<ncomp_to_dump_; i++) + { + HeaderFile << mesh_field_names_[i] << ' '; + } + HeaderFile << "\n"; } #endif } diff --git a/Source/Diagnostics/BoostedFrame_module.F90 b/Source/Diagnostics/BoostedFrame_module.F90 deleted file mode 100644 index 3bce1817e..000000000 --- a/Source/Diagnostics/BoostedFrame_module.F90 +++ /dev/null @@ -1,106 +0,0 @@ - -module warpx_boosted_frame_module - - use iso_c_binding - use amrex_fort_module, only : amrex_real - use constants, only : clight - - implicit none - -contains - -! -! Given cell-centered data in the boosted reference frame of the simulation, -! this transforms E and B in place so that the multifab now contains values -! in the lab frame. This routine assumes that the simulation frame is moving -! in the positive z direction with respect to the lab frame. -! - subroutine warpx_lorentz_transform_z(data, dlo, dhi, tlo, thi, gamma_boost, beta_boost) & - bind(C, name="warpx_lorentz_transform_z") - - integer(c_int), intent(in) :: dlo(3), dhi(3) - integer(c_int), intent(in) :: tlo(3), thi(3) - real(amrex_real), intent(inout) :: data(dlo(1):dhi(1),dlo(2):dhi(2),dlo(3):dhi(3), 10) - real(amrex_real), intent(in) :: gamma_boost, beta_boost - - integer i, j, k - real(amrex_real) e_lab, b_lab, j_lab, r_lab - - do k = tlo(3), thi(3) - do j = tlo(2), thi(2) - do i = tlo(1), thi(1) - - ! Transform the transverse E and B fields. Note that ez and bz are not - ! changed by the tranform. - e_lab = gamma_boost * (data(i, j, k, 1) + beta_boost*clight*data(i, j, k, 5)) - b_lab = gamma_boost * (data(i, j, k, 5) + beta_boost*data(i, j, k, 1)/clight) - - data(i, j, k, 1) = e_lab - data(i, j, k, 5) = b_lab - - e_lab = gamma_boost * (data(i, j, k, 2) - beta_boost*clight*data(i, j, k, 4)) - b_lab = gamma_boost * (data(i, j, k, 4) - beta_boost*data(i, j, k, 2)/clight) - - data(i, j, k, 2) = e_lab - data(i, j, k, 4) = b_lab - - ! Transform the charge and current density. Only the z component of j is affected. - j_lab = gamma_boost*(data(i, j, k, 9) + beta_boost*clight*data(i, j, k, 10)) - r_lab = gamma_boost*(data(i, j, k, 10) + beta_boost*data(i, j, k, 9)/clight) - - data(i, j, k, 9) = j_lab - data(i, j, k, 10) = r_lab - - end do - end do - end do - - end subroutine warpx_lorentz_transform_z - - subroutine warpx_copy_slice_3d(lo, hi, tmp, tlo, thi, & - buf, blo, bhi, ncomp, & - i_boost, i_lab) & - bind(C, name="warpx_copy_slice_3d") - - integer , intent(in) :: ncomp, i_boost, i_lab - integer , intent(in) :: lo(3), hi(3) - integer , intent(in) :: tlo(3), thi(3) - integer , intent(in) :: blo(3), bhi(3) - real(amrex_real), intent(inout) :: tmp(tlo(1):thi(1),tlo(2):thi(2),tlo(3):thi(3),ncomp) - real(amrex_real), intent(inout) :: buf(blo(1):bhi(1),blo(2):bhi(2),blo(3):bhi(3),ncomp) - - integer n, i, j, k - - do n = 1, ncomp - do j = lo(2), hi(2) - do i = lo(1), hi(1) - buf(i, j, i_lab, n) = tmp(i, j, i_boost, n) - end do - end do - end do - - end subroutine warpx_copy_slice_3d - - subroutine warpx_copy_slice_2d(lo, hi, tmp, tlo, thi, & - buf, blo, bhi, ncomp, & - i_boost, i_lab) & - bind(C, name="warpx_copy_slice_2d") - - integer , intent(in) :: ncomp, i_boost, i_lab - integer , intent(in) :: lo(2), hi(2) - integer , intent(in) :: tlo(2), thi(2) - integer , intent(in) :: blo(2), bhi(2) - real(amrex_real), intent(inout) :: tmp(tlo(1):thi(1),tlo(2):thi(2),ncomp) - real(amrex_real), intent(inout) :: buf(blo(1):bhi(1),blo(2):bhi(2),ncomp) - - integer n, i, j, k - - do n = 1, ncomp - do i = lo(1), hi(1) - buf(i, i_lab, n) = tmp(i, i_boost, n) - end do - end do - - end subroutine warpx_copy_slice_2d - -end module warpx_boosted_frame_module diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index 209d8e9b4..b1181f22f 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -299,7 +299,10 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, amrex::Vector<MultiFab>& mf_avg, const int ngrow) const { // Count how many different fields should be written (ncomp) - const int ncomp = 3*3 + const int ncomp = 0 + + static_cast<int>(plot_E_field)*3 + + static_cast<int>(plot_B_field)*3 + + static_cast<int>(plot_J_field)*3 + static_cast<int>(plot_part_per_cell) + static_cast<int>(plot_part_per_grid) + static_cast<int>(plot_part_per_proc) @@ -321,15 +324,21 @@ WarpX::AverageAndPackFields ( Vector<std::string>& varnames, // Go through the different fields, pack them into mf_avg[lev], // add the corresponding names to `varnames` and increment dcomp int dcomp = 0; - AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"jx","jy","jz"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"Ex","Ey","Ez"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"Bx","By","Bz"}) varnames.push_back(name); - dcomp += 3; + if (plot_J_field) { + AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"jx","jy","jz"}) varnames.push_back(name); + dcomp += 3; + } + if (plot_E_field) { + AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"Ex","Ey","Ez"}) varnames.push_back(name); + dcomp += 3; + } + if (plot_B_field) { + AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"Bx","By","Bz"}) varnames.push_back(name); + dcomp += 3; + } if (plot_part_per_cell) { @@ -651,7 +660,6 @@ getInterpolatedScalar( interpolated_F->setVal(0.); // Loop through the boxes and interpolate the values from the _cp data - const int use_limiter = 0; #ifdef _OPEMP #pragma omp parallel #endif @@ -669,7 +677,7 @@ getInterpolatedScalar( if ( F_fp.is_nodal() ){ IntVect refinement_vector{AMREX_D_DECL(r_ratio, r_ratio, r_ratio)}; node_bilinear_interp.interp(cfab, 0, ffab, 0, 1, - finebx, refinement_vector, {}, {}, {}, 0, 0); + finebx, refinement_vector, {}, {}, {}, 0, 0, RunOn::Cpu); } else { amrex::Abort("Unknown field staggering."); } diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package index 1c602ee33..dfd947d53 100644 --- a/Source/Diagnostics/Make.package +++ b/Source/Diagnostics/Make.package @@ -5,7 +5,8 @@ CEXE_sources += FieldIO.cpp CEXE_headers += FieldIO.H CEXE_headers += BoostedFrameDiagnostic.H CEXE_headers += ElectrostaticIO.cpp -F90EXE_sources += BoostedFrame_module.F90 +CEXE_headers += SliceDiagnostic.H +CEXE_sources += SliceDiagnostic.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index f15c084a0..f2a543ed5 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -5,6 +5,53 @@ using namespace amrex; void +RigidInjectedParticleContainer::ReadHeader (std::istream& is) +{ + is >> charge >> mass; + WarpX::GotoNextLine(is); + + int nlevs; + is >> nlevs; + WarpX::GotoNextLine(is); + + AMREX_ASSERT(zinject_plane_levels.size() == 0); + AMREX_ASSERT(done_injecting.size() == 0); + + for (int i = 0; i < nlevs; ++i) + { + int zinject_plane_tmp; + is >> zinject_plane_tmp; + zinject_plane_levels.push_back(zinject_plane_tmp); + WarpX::GotoNextLine(is); + } + + for (int i = 0; i < nlevs; ++i) + { + int done_injecting_tmp; + is >> done_injecting_tmp; + done_injecting.push_back(done_injecting_tmp); + WarpX::GotoNextLine(is); + } +} + +void +RigidInjectedParticleContainer::WriteHeader (std::ostream& os) const +{ + // no need to write species_id + os << charge << " " << mass << "\n"; + int nlevs = zinject_plane_levels.size(); + os << nlevs << "\n"; + for (int i = 0; i < nlevs; ++i) + { + os << zinject_plane_levels[i] << "\n"; + } + for (int i = 0; i < nlevs; ++i) + { + os << done_injecting[i] << "\n"; + } +} + +void WarpXParticleContainer::ReadHeader (std::istream& is) { is >> charge >> mass; @@ -27,17 +74,55 @@ MultiParticleContainer::Checkpoint (const std::string& dir) const } void -MultiParticleContainer::WritePlotFile (const std::string& dir, - const Vector<int>& real_flags, - const Vector<std::string>& real_names) const +MultiParticleContainer::WritePlotFile (const std::string& dir) const { Vector<std::string> int_names; Vector<int> int_flags; - + for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - allcontainers[i]->WritePlotFile(dir, species_names[i], - real_flags, int_flags, - real_names, int_names); + auto& pc = allcontainers[i]; + if (pc->plot_species) { + + Vector<std::string> real_names; + real_names.push_back("weight"); + + 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_RZ + real_names.push_back("theta"); +#endif + + if (WarpX::do_boosted_frame_diagnostic && pc->DoBoostedFrameDiags()) + { + real_names.push_back("xold"); + real_names.push_back("yold"); + real_names.push_back("zold"); + + real_names.push_back("uxold"); + real_names.push_back("uyold"); + real_names.push_back("uzold"); + } + + // Convert momentum to SI + pc->ConvertUnits(ConvertDirection::WarpX_to_SI); + // real_names contains a list of all particle attributes. + // pc->plot_flags is 1 or 0, whether quantity is dumped or not. + pc->WritePlotFile(dir, species_names[i], + pc->plot_flags, int_flags, + real_names, int_names); + // Convert momentum back to WarpX units + pc->ConvertUnits(ConvertDirection::SI_to_WarpX); + } } } @@ -65,3 +150,43 @@ MultiParticleContainer::WriteHeader (std::ostream& os) const } } +// 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. +void +PhysicalParticleContainer::ConvertUnits(ConvertDirection convert_direction) +{ + BL_PROFILE("PPC::ConvertUnits()"); + + // Compute conversion factor + Real factor = 1; + if (convert_direction == ConvertDirection::WarpX_to_SI){ + factor = mass; + } else if (convert_direction == ConvertDirection::SI_to_WarpX){ + factor = 1./mass; + } + + for (int lev=0; lev<=finestLevel(); lev++){ +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { + // - momenta are stored as a struct of array, in `attribs` + auto& attribs = pti.GetAttribs(); + Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + // Loop over the particles and convert momentum + const long np = pti.numParticles(); + ParallelFor( np, + [=] AMREX_GPU_DEVICE (long i) { + ux[i] *= factor; + uy[i] *= factor; + uz[i] *= factor; + } + ); + } + } +} diff --git a/Source/Diagnostics/SliceDiagnostic.H b/Source/Diagnostics/SliceDiagnostic.H new file mode 100644 index 000000000..31eea83be --- /dev/null +++ b/Source/Diagnostics/SliceDiagnostic.H @@ -0,0 +1,42 @@ +#ifndef WARPX_SliceDiagnostic_H_ +#define WARPX_SliceDiagnostic_H_ + +#include <AMReX_VisMF.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_ParallelDescriptor.H> +#include <AMReX_Geometry.H> + +#include <WarpX.H> +#include <AMReX_FArrayBox.H> +#include <AMReX_IArrayBox.H> +#include <AMReX_Vector.H> +#include <AMReX_BLassert.H> +#include <AMReX_MultiFabUtil.H> +#include <AMReX_MultiFabUtil_C.H> + + + +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, + amrex::IntVect &slice_cr_ratio, amrex::Vector<amrex::Geometry> dom_geom, + 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, + 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, + int srccomp, int ncomp, int nghost, const amrex::RealBox real_box); + +#endif + diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp new file mode 100644 index 000000000..994f990c6 --- /dev/null +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -0,0 +1,475 @@ +#include "SliceDiagnostic.H" +#include <AMReX_MultiFabUtil.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_FillPatchUtil_F.H> + +#include <WarpX.H> + +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 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 + */ + +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; + std::unique_ptr<MultiFab> cs_mf; + + Vector<int> slice_ncells(AMREX_SPACEDIM); + int nghost = 1; + int nlevels = dom_geom.size(); + int ncomp = (mf).nComp(); + + 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 ) + { + SliceType[idim] = conversionType.nodeCentered(idim); + } + + const RealBox& real_box = dom_geom[0].ProbDomain(); + RealBox slice_cc_nd_box; + int slice_grid_size = 32; + + bool interpolate = false; + bool coarsen = false; + + // same index space as domain // + IntVect slice_lo(AMREX_D_DECL(0,0,0)); + IntVect slice_hi(AMREX_D_DECL(1,1,1)); + 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, + slice_hi, interp_lo); + // Determine if interpolation is required and number of cells in slice // + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + + // Flag for interpolation if required // + if ( interp_lo[idim] == 1) { + interpolate = 1; + } + + // For the case when a dimension is reduced // + if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { + slice_ncells[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) { + coarsen = true; + + // modify slice_grid_size if >= refines_cells // + if ( slice_grid_size >= refined_ncells ) { + slice_grid_size = refined_ncells - 1; + } + + } + } + } + + // 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); + + // 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)); + + // Copy data from domain to slice that has same cell size as that of // + // the domain mf. src and dst have the same number of ghost cells // + 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, + ncomp, nghost, slice_lo, slice_hi, SliceType, real_box); + } + + + if (coarsen == false) { + return smf; + } + else if ( coarsen == true ) { + Vector<BoxArray> crse_ba(1); + crse_ba[0] = sba[0]; + crse_ba[0].coarsen(slice_cr_ratio); + + AMREX_ALWAYS_ASSERT(crse_ba[0].size() == sba[0].size()); + + 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) { + + FArrayBox& Src_fabox = mfSrc[mfi]; + + const Box& Dst_bx = mfi_dst.validbox(); + FArrayBox& Dst_fabox = mfDst[mfi_dst]; + + int scomp = 0; + int dcomp = 0; + + IntVect cctype(AMREX_D_DECL(0,0,0)); + if( SliceType==cctype ) { + amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, + ncomp, slice_cr_ratio); + } + IntVect ndtype(AMREX_D_DECL(1,1,1)); + if( SliceType == ndtype ) { + amrex::amrex_avgdown_nodes(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio); + } + if( SliceType == WarpX::Ex_nodal_flag ) { + 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, + scomp, ncomp, slice_cr_ratio, 1); + } + if( SliceType == WarpX::Ez_nodal_flag ) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 2); + } + if( SliceType == WarpX::Bx_nodal_flag) { + 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, + scomp, ncomp, slice_cr_ratio, 1); + } + if( SliceType == WarpX::Bz_nodal_flag ) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 2); + } + + if ( mfi_dst.isValid() ) { + ++mfi_dst; + } + + } + return cs_mf; + + } + amrex::Abort("Should not hit this return statement."); + return smf; +} + + +/* \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 + * 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 + * 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. + * \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. + */ + + +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) + { + // Modify coarsening ratio if the input value is not an exponent of 2 for AMR // + 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) ) { + slice_realbox.setLo( idim, real_box.lo(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) ) { + slice_realbox.setHi( idim, real_box.hi(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 // + if ( slice_realbox.hi(idim) - slice_realbox.lo(idim) <= 0) + { + 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 ) ) + / 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) ) ) + / dom_geom[0].CellSize(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_cc_nd_box.lo(idim) - real_box.lo(idim) < fac ) { + interp_lo[idim] = 1; + } + } + 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 + - (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) { + int lo_new = index_lo; + 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; + } + // To ensure that the index.hi is coarsenable // + if ( mod_hi > 0) { + hi_new = index_hi + (slice_cr_ratio[idim] - mod_hi); + } + + // If modified index.hi is > baselinebox.hi, move the point // + // to the previous coarsenable point // + 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 ){ + 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 ) { + 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_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) + + 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 ); + } + } +} + + +/* \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. + */ +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, + const RealBox real_box) +{ + for (MFIter mfi(smf); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const auto IndType = smf.ixType(); + const auto lo = amrex::lbound(bx); + 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, + slice_realbox, 0, ncomp, nghost, real_box); + } + } + } + +} + +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(); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + double fac = ( 1.0-IndType[idir] )*geom[0].CellSize(idir) * 0.5; + int imin = slice_lo[idir]; + double minpos = imin*geom[0].CellSize(idir) + fac + real_box.lo(idir); + double maxpos = (imin+1)*geom[0].CellSize(idir) + fac + real_box.lo(idir); + double slice_minpos = slice_realbox.lo(idir) ; + + switch (idir) { + case 0: + { + if ( imin >= lo.x && imin <= lo.x) { + for (int n = srccomp; n < srccomp + ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i+1,j,k,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + 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) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i,j+1,k,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + case 2: + { + if ( imin >= lo.z && imin <= lo.z) { + for (int n = srccomp; n < srccomp+ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i,j,k+1,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + + } + +} + + + + + + + diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 186a8d45e..869d3580e 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -11,6 +11,8 @@ #include <AMReX_AmrMeshInSituBridge.H> #endif +#include "SliceDiagnostic.H" + #ifdef AMREX_USE_ASCENT #include <ascent.hpp> #include <AMReX_Conduit_Blueprint.H> @@ -84,11 +86,11 @@ WarpX::WriteWarpXHeader(const std::string& name) const // Geometry for (int i = 0; i < AMREX_SPACEDIM; ++i) { - HeaderFile << Geometry::ProbLo(i) << ' '; + HeaderFile << Geom(0).ProbLo(i) << ' '; } HeaderFile << '\n'; for (int i = 0; i < AMREX_SPACEDIM; ++i) { - HeaderFile << Geometry::ProbHi(i) << ' '; + HeaderFile << Geom(0).ProbHi(i) << ' '; } HeaderFile << '\n'; @@ -283,7 +285,7 @@ WarpX::InitFromCheckpoint () } } - Geometry::ProbDomain(RealBox(prob_lo,prob_hi)); + ResetProbDomain(RealBox(prob_lo,prob_hi)); for (int lev = 0; lev < nlevs; ++lev) { BoxArray ba; @@ -422,14 +424,14 @@ WarpX::GetCellCenteredData() { AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dcomp, ng ); dcomp += 3; // then the magnetic field - AverageAndPackVectorField( *cc[lev], Efield_aux[lev], dcomp, ng ); + AverageAndPackVectorField( *cc[lev], Bfield_aux[lev], dcomp, ng ); dcomp += 3; // then the current density AverageAndPackVectorField( *cc[lev], current_fp[lev], dcomp, ng ); dcomp += 3; + // then the charge density const std::unique_ptr<MultiFab>& charge_density = mypc->GetChargeDensity(lev); AverageAndPackScalarField( *cc[lev], *charge_density, dcomp, ng ); - cc[lev]->FillBoundary(geom[lev].periodicity()); } @@ -612,37 +614,7 @@ WarpX::WritePlotFile () const } } - Vector<std::string> particle_varnames; - particle_varnames.push_back("weight"); - - particle_varnames.push_back("momentum_x"); - particle_varnames.push_back("momentum_y"); - particle_varnames.push_back("momentum_z"); - - particle_varnames.push_back("Ex"); - particle_varnames.push_back("Ey"); - particle_varnames.push_back("Ez"); - - particle_varnames.push_back("Bx"); - particle_varnames.push_back("By"); - particle_varnames.push_back("Bz"); - -#ifdef WARPX_RZ - particle_varnames.push_back("theta"); -#endif - - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) - { - particle_varnames.push_back("xold"); - particle_varnames.push_back("yold"); - particle_varnames.push_back("zold"); - - particle_varnames.push_back("uxold"); - particle_varnames.push_back("uyold"); - particle_varnames.push_back("uzold"); - } - - mypc->WritePlotFile(plotfilename, particle_plot_flags, particle_varnames); + mypc->WritePlotFile(plotfilename); WriteJobInfo(plotfilename); @@ -771,3 +743,94 @@ WarpX::WriteJobInfo (const std::string& dir) const jobInfoFile.close(); } } + + +/* \brief + * The slice is ouput using visMF and can be visualized used amrvis. + */ +void +WarpX::WriteSlicePlotFile () const +{ + if (F_fp[0] ) { + VisMF::Write( (*F_slice[0]), "vismf_F_slice"); + } + + if (rho_fp[0]) { + VisMF::Write( (*rho_slice[0]), "vismf_rho_slice"); + } + + VisMF::Write( (*Efield_slice[0][0]), amrex::Concatenate("vismf_Ex_slice_",istep[0])); + VisMF::Write( (*Efield_slice[0][1]), amrex::Concatenate("vismf_Ey_slice_",istep[0])); + VisMF::Write( (*Efield_slice[0][2]), amrex::Concatenate("vismf_Ez_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][0]), amrex::Concatenate("vismf_Bx_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][1]), amrex::Concatenate("vismf_By_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][2]), amrex::Concatenate("vismf_Bz_slice_",istep[0])); + VisMF::Write( (*current_slice[0][0]), amrex::Concatenate("vismf_jx_slice_",istep[0])); + VisMF::Write( (*current_slice[0][1]), amrex::Concatenate("vismf_jy_slice_",istep[0])); + VisMF::Write( (*current_slice[0][2]), amrex::Concatenate("vismf_jz_slice_",istep[0])); + +} + + +void +WarpX::InitializeSliceMultiFabs () +{ + + int nlevels = Geom().size(); + + F_slice.resize(nlevels); + rho_slice.resize(nlevels); + current_slice.resize(nlevels); + Efield_slice.resize(nlevels); + Bfield_slice.resize(nlevels); + +} + + +// To generate slice that inherits index type of underlying data // +void +WarpX::SliceGenerationForDiagnostics () +{ + + Vector<Geometry> dom_geom; + dom_geom = Geom(); + + if (F_fp[0] ) { + 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, + slice_cr_ratio ); + } + + for (int idim = 0; idim < 3; ++idim) { + 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(), + dom_geom, slice_realbox, slice_cr_ratio ); + current_slice[0][idim] = CreateSlice( *current_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + } + + +} + + +void +WarpX::ClearSliceMultiFabs () +{ + + F_slice.clear(); + rho_slice.clear(); + current_slice.clear(); + Efield_slice.clear(); + Bfield_slice.clear(); + F_slice.shrink_to_fit(); + rho_slice.shrink_to_fit(); + current_slice.shrink_to_fit(); + Efield_slice.shrink_to_fit(); + Bfield_slice.shrink_to_fit(); + +} + diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index e98561be1..4f33694cd 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -25,6 +25,10 @@ WarpX::EvolveEM (int numsteps) static int last_check_file_step = 0; static int last_insitu_step = 0; + if (do_compute_max_step_from_zmax){ + computeMaxStepBoostAccelerator(geom[0]); + } + int numsteps_max; if (numsteps < 0) { // Note that the default argument is numsteps = -1 numsteps_max = max_step; @@ -80,12 +84,14 @@ WarpX::EvolveEM (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1],*Bfield_aux[lev][2]); } is_synchronized = false; + } else { // Beyond one step, we have E^{n} and B^{n}. // Particles have p^{n-1/2} and x^{n}. FillBoundaryE(); FillBoundaryB(); UpdateAuxilaryData(); + } if (do_subcycling == 0 || finest_level == 0) { @@ -128,6 +134,8 @@ WarpX::EvolveEM (int numsteps) bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); + // slice generation // + bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); @@ -172,7 +180,8 @@ WarpX::EvolveEM (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - if (to_make_plot || do_insitu) + // slice gen // + if (to_make_plot || do_insitu || to_make_slice_plot) { FillBoundaryE(); FillBoundaryB(); @@ -188,11 +197,19 @@ WarpX::EvolveEM (int numsteps) last_insitu_step = step+1; if (to_make_plot) - WritePlotFile(); + WritePlotFile(); + + if (to_make_slice_plot) + { + InitializeSliceMultiFabs (); + SliceGenerationForDiagnostics(); + WriteSlicePlotFile(); + ClearSliceMultiFabs (); + } if (do_insitu) UpdateInSitu(); - } + } if (check_int > 0 && (step+1) % check_int == 0) { last_check_file_step = step+1; @@ -268,6 +285,7 @@ WarpX::OneStep_nosub (Real cur_time) if (warpx_py_beforedeposition) warpx_py_beforedeposition(); #endif PushParticlesandDepose(cur_time); + #ifdef WARPX_USE_PY if (warpx_py_afterdeposition) warpx_py_afterdeposition(); #endif @@ -503,6 +521,50 @@ WarpX::ComputeDt () } } +/* \brief computes max_step for wakefield simulation in boosted frame. + * \param geom: Geometry object that contains simulation domain. + * + * max_step is set so that the simulation stop when the lower corner of the + * simulation box passes input parameter zmax_plasma_to_compute_max_step. + */ +void +WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ + // Sanity checks: can use zmax_plasma_to_compute_max_step only if + // the moving window and the boost are all in z direction. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + WarpX::moving_window_dir == AMREX_SPACEDIM-1, + "Can use zmax_plasma_to_compute_max_step only if " + + "moving window along z. TODO: all directions."); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + maxLevel() == 0, + "Can use zmax_plasma_to_compute_max_step only if " + + "max level = 0."); + 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[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, + "Can use zmax_plasma_to_compute_max_step only if " + + "warpx.boost_direction = z. TODO: all directions."); + + // Lower end of the simulation domain. All quantities are given in boosted + // frame except zmax_plasma_to_compute_max_step. + const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1); + // End of the plasma: Transform input argument + // zmax_plasma_to_compute_max_step to boosted frame. + const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; + // Plasma velocity + const Real v_plasma_boost = -beta_boost * PhysConst::c; + // Get time at which the lower end of the simulation domain passes the + // upper end of the plasma (in the z direction). + const Real interaction_time_boost = (len_plasma_boost-zmin_domain_boost)/ + (moving_window_v-v_plasma_boost); + // Divide by dt, and update value of max_step. + const int computed_max_step = interaction_time_boost/dt[0]; + max_step = computed_max_step; + Print()<<"max_step computed in computeMaxStepBoostAccelerator: " + <<computed_max_step<<std::endl; +} + /* \brief Apply perfect mirror condition inside the box (not at a boundary). * In practice, set all fields to 0 on a section of the simulation domain * (as for a perfect conductor with a given thickness). @@ -543,19 +605,19 @@ WarpX::applyMirrors(Real time){ NullifyMF(Bz, lev, z_min, z_max); if (lev>0){ // Get coarse patch field MultiFabs - MultiFab& Ex = *Efield_cp[lev][0].get(); - MultiFab& Ey = *Efield_cp[lev][1].get(); - MultiFab& Ez = *Efield_cp[lev][2].get(); - MultiFab& Bx = *Bfield_cp[lev][0].get(); - MultiFab& By = *Bfield_cp[lev][1].get(); - MultiFab& Bz = *Bfield_cp[lev][2].get(); + MultiFab& cEx = *Efield_cp[lev][0].get(); + MultiFab& cEy = *Efield_cp[lev][1].get(); + MultiFab& cEz = *Efield_cp[lev][2].get(); + MultiFab& cBx = *Bfield_cp[lev][0].get(); + MultiFab& cBy = *Bfield_cp[lev][1].get(); + MultiFab& cBz = *Bfield_cp[lev][2].get(); // Set each field to zero between z_min and z_max - NullifyMF(Ex, lev, z_min, z_max); - NullifyMF(Ey, lev, z_min, z_max); - NullifyMF(Ez, lev, z_min, z_max); - NullifyMF(Bx, lev, z_min, z_max); - NullifyMF(By, lev, z_min, z_max); - NullifyMF(Bz, lev, z_min, z_max); + NullifyMF(cEx, lev, z_min, z_max); + NullifyMF(cEy, lev, z_min, z_max); + NullifyMF(cEz, lev, z_min, z_max); + NullifyMF(cBx, lev, z_min, z_max); + NullifyMF(cBy, lev, z_min, z_max); + NullifyMF(cBz, lev, z_min, z_max); } } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H index 0487e5226..12718e38b 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H @@ -8,14 +8,18 @@ */ class PsatdAlgorithm : public SpectralBaseAlgorithm { + public: PsatdAlgorithm(const SpectralKSpace& spectral_kspace, const amrex::DistributionMapping& dm, const int norder_x, const int norder_y, - const int norder_z, const bool nodal, - const amrex::Real dt); - // Redefine update equation from base class - virtual void pushSpectralFields(SpectralFieldData& f) const override final; + const int norder_z, const bool nodal, const amrex::Real dt); + + void InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt); + + void pushSpectralFields(SpectralFieldData& f) const override final; private: SpectralCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp index 37892d35a..d45b01bda 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp @@ -22,59 +22,8 @@ PsatdAlgorithm::PsatdAlgorithm(const SpectralKSpace& spectral_kspace, X2_coef = SpectralCoefficients(ba, dm, 1, 0); X3_coef = SpectralCoefficients(ba, dm, 1, 0); - // Fill them with the right values: - // Loop over boxes and allocate the corresponding coefficients - // for each box owned by the local MPI proc - for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ - - const Box& bx = ba[mfi]; - - // Extract pointers for the k vectors - const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); -#if (AMREX_SPACEDIM==3) - const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); -#endif - const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); - // Extract arrays for the coefficients - Array4<Real> C = C_coef[mfi].array(); - Array4<Real> S_ck = S_ck_coef[mfi].array(); - Array4<Real> X1 = X1_coef[mfi].array(); - Array4<Real> X2 = X2_coef[mfi].array(); - Array4<Real> X3 = X3_coef[mfi].array(); - - // Loop over indices within one box - ParallelFor(bx, - [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept - { - // Calculate norm of vector - const Real k_norm = std::sqrt( - std::pow(modified_kx[i], 2) + -#if (AMREX_SPACEDIM==3) - std::pow(modified_ky[j], 2) + - std::pow(modified_kz[k], 2)); -#else - std::pow(modified_kz[j], 2)); -#endif - - // Calculate coefficients - constexpr Real c = PhysConst::c; - constexpr Real ep0 = PhysConst::ep0; - if (k_norm != 0){ - C(i,j,k) = std::cos(c*k_norm*dt); - S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); - X1(i,j,k) = (1. - C(i,j,k))/(ep0 * c*c * k_norm*k_norm); - X2(i,j,k) = (1. - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); - X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); - } else { // Handle k_norm = 0, by using the analytical limit - C(i,j,k) = 1.; - S_ck(i,j,k) = dt; - X1(i,j,k) = 0.5 * dt*dt / ep0; - X2(i,j,k) = c*c * dt*dt / (6.*ep0); - X3(i,j,k) = - c*c * dt*dt / (3.*ep0); - } - }); - } -}; + InitializeSpectralCoefficients(spectral_kspace, dm, dt); +} /* Advance the E and B field in spectral space (stored in `f`) * over one time step */ @@ -130,13 +79,14 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ #endif constexpr Real c2 = PhysConst::c*PhysConst::c; constexpr Real inv_ep0 = 1./PhysConst::ep0; - constexpr Complex I = Complex{0,1}; + const Complex I = Complex{0,1}; const Real C = C_arr(i,j,k); const Real S_ck = S_ck_arr(i,j,k); const Real X1 = X1_arr(i,j,k); const Real X2 = X2_arr(i,j,k); const Real X3 = X3_arr(i,j,k); + // Update E (see WarpX online documentation: theory section) fields(i,j,k,Idx::Ex) = C*Ex_old + S_ck*(c2*I*(ky*Bz_old - kz*By_old) - inv_ep0*Jx) @@ -160,3 +110,63 @@ PsatdAlgorithm::pushSpectralFields(SpectralFieldData& f) const{ }); } }; + +void PsatdAlgorithm::InitializeSpectralCoefficients(const SpectralKSpace& spectral_kspace, + const amrex::DistributionMapping& dm, + const amrex::Real dt) +{ + const BoxArray& ba = spectral_kspace.spectralspace_ba; + // Fill them with the right values: + // Loop over boxes and allocate the corresponding coefficients + // for each box owned by the local MPI proc + for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi){ + + const Box& bx = ba[mfi]; + + // Extract pointers for the k vectors + const Real* modified_kx = modified_kx_vec[mfi].dataPtr(); +#if (AMREX_SPACEDIM==3) + const Real* modified_ky = modified_ky_vec[mfi].dataPtr(); +#endif + const Real* modified_kz = modified_kz_vec[mfi].dataPtr(); + // Extract arrays for the coefficients + Array4<Real> C = C_coef[mfi].array(); + Array4<Real> S_ck = S_ck_coef[mfi].array(); + Array4<Real> X1 = X1_coef[mfi].array(); + Array4<Real> X2 = X2_coef[mfi].array(); + Array4<Real> X3 = X3_coef[mfi].array(); + + // Loop over indices within one box + ParallelFor(bx, + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept + { + // Calculate norm of vector + const Real k_norm = std::sqrt( + std::pow(modified_kx[i], 2) + +#if (AMREX_SPACEDIM==3) + std::pow(modified_ky[j], 2) + + std::pow(modified_kz[k], 2)); +#else + std::pow(modified_kz[j], 2)); +#endif + + + // Calculate coefficients + constexpr Real c = PhysConst::c; + constexpr Real ep0 = PhysConst::ep0; + if (k_norm != 0){ + C(i,j,k) = std::cos(c*k_norm*dt); + S_ck(i,j,k) = std::sin(c*k_norm*dt)/(c*k_norm); + X1(i,j,k) = (1. - C(i,j,k))/(ep0 * c*c * k_norm*k_norm); + X2(i,j,k) = (1. - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); + X3(i,j,k) = (C(i,j,k) - S_ck(i,j,k)/dt)/(ep0 * k_norm*k_norm); + } else { // Handle k_norm = 0, by using the analytical limit + C(i,j,k) = 1.; + S_ck(i,j,k) = dt; + X1(i,j,k) = 0.5 * dt*dt / ep0; + X2(i,j,k) = c*c * dt*dt / (6.*ep0); + X3(i,j,k) = - c*c * dt*dt / (3.*ep0); + } + }); + } +} diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 8e58aa1d8..7954414b8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -25,7 +25,7 @@ class SpectralFieldData // (plans are only initialized for the boxes that are owned by // the local MPI rank) #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + using FFTplans = amrex::LayoutData<cufftHandle>; #else using FFTplans = amrex::LayoutData<fftw_plan>; #endif diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 02fa2015f..a2b695568 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -53,7 +53,38 @@ SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, // the FFT plan, the valid dimensions are those of the real-space box. IntVect fft_size = realspace_ba[mfi].length(); #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + // Create cuFFT plans + // Creating 3D plan for real to complex -- double precision + // Assuming CUDA is used for programming GPU + // Note that D2Z is inherently forward plan + // and Z2D is inherently backward plan + cufftResult result; +#if (AMREX_SPACEDIM == 3) + result = cufftPlan3d( &forward_plan[mfi], fft_size[2], + fft_size[1],fft_size[0], CUFFT_D2Z); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " cufftplan3d forward failed! \n"; + } + + result = cufftPlan3d( &backward_plan[mfi], fft_size[2], + fft_size[1], fft_size[0], CUFFT_Z2D); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " cufftplan3d backward failed! \n"; + } +#else + result = cufftPlan2d( &forward_plan[mfi], fft_size[1], + fft_size[0], CUFFT_D2Z ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " cufftplan2d forward failed! \n"; + } + + result = cufftPlan2d( &backward_plan[mfi], fft_size[1], + fft_size[0], CUFFT_Z2D ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " cufftplan2d backward failed! \n"; + } +#endif + #else // Create FFTW plans forward_plan[mfi] = @@ -86,7 +117,9 @@ SpectralFieldData::~SpectralFieldData() if (tmpRealField.size() > 0){ for ( MFIter mfi(tmpRealField); mfi.isValid(); ++mfi ){ #ifdef AMREX_USE_GPU - // Add cuFFT-specific code + // Destroy cuFFT plans + cufftDestroy( forward_plan[mfi] ); + cufftDestroy( backward_plan[mfi] ); #else // Destroy FFTW plans fftw_destroy_plan( forward_plan[mfi] ); @@ -129,14 +162,25 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, Array4<Real> tmp_arr = tmpRealField[mfi].array(); ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { - tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); + tmp_arr(i,j,k) = mf_arr(i,j,k,i_comp); }); } // Perform Fourier transform from `tmpRealField` to `tmpSpectralField` #ifdef AMREX_USE_GPU - // Add cuFFT-specific code ; make sure that this is done on the same - // GPU stream as the above copy + // Perform Fast Fourier Transform on GPU using cuFFT + // make sure that this is done on the same + // GPU stream as the above copy + cufftResult result; + cudaStream_t stream = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( forward_plan[mfi], stream); + result = cufftExecD2Z( forward_plan[mfi], + tmpRealField[mfi].dataPtr(), + reinterpret_cast<cuDoubleComplex*>( + tmpSpectralField[mfi].dataPtr()) ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " forward transform using cufftExecD2Z failed ! \n"; + } #else fftw_execute( forward_plan[mfi] ); #endif @@ -155,6 +199,7 @@ SpectralFieldData::ForwardTransform( const MultiFab& mf, const Complex* zshift_arr = zshift_FFTfromCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = tmp_arr(i,j,k); @@ -207,6 +252,7 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, const Complex* zshift_arr = zshift_FFTtoCell[mfi].dataPtr(); // Loop over indices within one box const Box spectralspace_bx = tmpSpectralField[mfi].box(); + ParallelFor( spectralspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { Complex spectral_field_value = field_arr(i,j,k,field_index); @@ -225,8 +271,19 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, // Perform Fourier transform from `tmpSpectralField` to `tmpRealField` #ifdef AMREX_USE_GPU - // Add cuFFT-specific code ; make sure that this is done on the same + // Perform Fast Fourier Transform on GPU using cuFFT. + // make sure that this is done on the same // GPU stream as the above copy + cufftResult result; + cudaStream_t stream = amrex::Gpu::Device::cudaStream(); + cufftSetStream ( backward_plan[mfi], stream); + result = cufftExecZ2D( backward_plan[mfi], + reinterpret_cast<cuDoubleComplex*>( + tmpSpectralField[mfi].dataPtr()), + tmpRealField[mfi].dataPtr() ); + if ( result != CUFFT_SUCCESS ) { + amrex::Print() << " Backward transform using cufftexecZ2D failed! \n"; + } #else fftw_execute( backward_plan[mfi] ); #endif @@ -240,6 +297,7 @@ SpectralFieldData::BackwardTransform( MultiFab& mf, Array4<const Real> tmp_arr = tmpRealField[mfi].array(); // Normalization: divide by the number of points in realspace const Real inv_N = 1./realspace_bx.numPts(); + ParallelFor( realspace_bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { // Copy and normalize field diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 2fe78cedd..6fe5e3939 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -142,9 +142,14 @@ SpectralKSpace::getSpectralShiftFactor( const DistributionMapping& dm, case ShiftType::TransformFromCellCentered: sign = -1.; break; case ShiftType::TransformToCellCentered: sign = 1.; } - constexpr Complex I{0,1}; + const Complex I{0,1}; for (int i=0; i<k.size(); i++ ){ +#ifdef AMREX_USE_GPU + shift[i] = thrust::exp( I*sign*k[i]*0.5*dx[i_dim] ); +#else shift[i] = std::exp( I*sign*k[i]*0.5*dx[i_dim] ); +#endif + } } return shift_factor; diff --git a/Source/FieldSolver/WarpXFFT.cpp b/Source/FieldSolver/WarpXFFT.cpp index 2047a569d..9815d43dc 100644 --- a/Source/FieldSolver/WarpXFFT.cpp +++ b/Source/FieldSolver/WarpXFFT.cpp @@ -56,6 +56,7 @@ BuildFFTOwnerMask (const MultiFab& mf, const Geometry& geom) for (const auto& b : bl) { fab.setVal(nonowner, b, 0, 1); } + } return mask; @@ -107,6 +108,8 @@ CopyDataFromFFTToValid (MultiFab& mf, const MultiFab& mf_fft, const BoxArray& ba // the cell that has non-zero mask is the one which is retained. mf.setVal(0.0, 0); mf.ParallelAdd(mftmp); + + } } @@ -481,4 +484,5 @@ WarpX::PushPSATD_hybridFFT (int lev, amrex::Real /* dt */) { amrex::Abort("WarpX::PushPSATD: TODO"); } + } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6a9ceae5a..c53e13f8f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -17,30 +17,30 @@ using namespace amrex; void -WarpX::EvolveB (Real dt) +WarpX::EvolveB (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveB(lev, dt); + EvolveB(lev, a_dt); } } void -WarpX::EvolveB (int lev, Real dt) +WarpX::EvolveB (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); - EvolveB(lev, PatchType::fine, dt); + EvolveB(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveB(lev, PatchType::coarse, dt); + EvolveB(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array<Real,3>& dx = WarpX::CellSize(patch_level); - Real dtsdx = dt/dx[0], dtsdy = dt/dx[1], dtsdz = dt/dx[2]; + Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -65,6 +65,10 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) MultiFab* cost = costs[lev].get(); const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + // xmin is only used by the picsar kernel with cylindrical geometry, + // in which case it is actually rmin. + const Real xmin = Geom(0).ProbLo(0); + // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -77,10 +81,6 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) const Box& tby = mfi.tilebox(By_nodal_flag); const Box& tbz = mfi.tilebox(Bz_nodal_flag); - // xmin is only used by the picsar kernel with cylindrical geometry, - // in which case it is actually rmin. - const Real xmin = mfi.tilebox().smallEnd(0)*dx[0]; - if (do_nodal) { auto const& Bxfab = Bx->array(mfi); auto const& Byfab = By->array(mfi); @@ -164,30 +164,30 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveE (Real dt) +WarpX::EvolveE (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveE(lev, dt); + EvolveE(lev, a_dt); } } void -WarpX::EvolveE (int lev, Real dt) +WarpX::EvolveE (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); - EvolveE(lev, PatchType::fine, dt); + EvolveE(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveE(lev, PatchType::coarse, dt); + EvolveE(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { - const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; - const Real c2dt = (PhysConst::c*PhysConst::c) * dt; + const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; + const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array<Real,3>& dx = WarpX::CellSize(patch_level); @@ -224,6 +224,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) MultiFab* cost = costs[lev].get(); const IntVect& rr = (lev > 0) ? refRatio(lev-1) : IntVect::TheUnitVector(); + // xmin is only used by the picsar kernel with cylindrical geometry, + // in which case it is actually rmin. + const Real xmin = Geom(0).ProbLo(0); + // Loop through the grids, and over the tiles within each grid #ifdef _OPENMP #pragma omp parallel if (Gpu::notInLaunchRegion()) @@ -236,10 +240,6 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) const Box& tey = mfi.tilebox(Ey_nodal_flag); const Box& tez = mfi.tilebox(Ez_nodal_flag); - // xmin is only used by the picsar kernel with cylindrical geometry, - // in which case it is actually rmin. - const Real xmin = mfi.tilebox().smallEnd(0)*dx[0]; - if (do_nodal) { auto const& Exfab = Ex->array(mfi); auto const& Eyfab = Ey->array(mfi); @@ -358,27 +358,27 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveF (Real dt, DtType dt_type) +WarpX::EvolveF (Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt, dt_type); + EvolveF(lev, a_dt, a_dt_type); } } void -WarpX::EvolveF (int lev, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; - EvolveF(lev, PatchType::fine, dt, dt_type); - if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); + EvolveF(lev, PatchType::fine, a_dt, a_dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, a_dt, a_dt_type); } void -WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -388,7 +388,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); - const std::array<Real,3> dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + const std::array<Real,3> dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; MultiFab *Ex, *Ey, *Ez, *rho, *F; if (patch_type == PatchType::fine) @@ -408,12 +408,12 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) F = F_cp[lev].get(); } - const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; + const int rhocomp = (a_dt_type == DtType::FirstHalf) ? 0 : 1; MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + MultiFab::Saxpy(*F, a_dt, src, 0, 0, 1, 0); if (do_pml && pml[lev]->ok()) { diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index 9bade3c77..4fcfec14b 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -1,30 +1,17 @@ -#include <AMReX_MultiFab.H> -#include <AMReX_CudaContainers.H> +#include <Filter.H> -#ifndef WARPX_FILTER_H_ -#define WARPX_FILTER_H_ +#ifndef WARPX_BILIN_FILTER_H_ +#define WARPX_BILIN_FILTER_H_ -class BilinearFilter +class BilinearFilter : public Filter { public: - BilinearFilter () = default; + + BilinearFilter() = default; void ComputeStencils(); - void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); amrex::IntVect npass_each_dir; - amrex::IntVect stencil_length_each_dir; - - // public for cuda - void Filter(const amrex::Box& tbx, - amrex::Array4<amrex::Real const> const& tmp, - amrex::Array4<amrex::Real > const& dst, - int scomp, int dcomp, int ncomp); - -private: - - amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z; - amrex::Dim3 slen; }; -#endif +#endif // #ifndef WARPX_BILIN_FILTER_H_ diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index f6acaa5df..68ebde687 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -68,173 +68,3 @@ void BilinearFilter::ComputeStencils(){ slen.z = 1; #endif } - - -#ifdef AMREX_USE_CUDA - -void -BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - BL_PROFILE("BilinearFilter::ApplyStencil()"); - ncomp = std::min(ncomp, srcmf.nComp()); - - for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) - { - const auto& src = srcmf.array(mfi); - const auto& dst = dstmf.array(mfi); - const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - - // tmpfab has enough ghost cells for the stencil - FArrayBox tmp_fab(gbx,ncomp); - Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early - auto const& tmp = tmp_fab.array(); - - // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcmf[mfi].box(); - AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, - { - if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { - tmp(i,j,k,n) = src(i,j,k,n+scomp); - } else { - tmp(i,j,k,n) = 0.0; - } - }); - - // Apply filter - Filter(tbx, tmp, dst, 0, dcomp, ncomp); - } -} - -void BilinearFilter::Filter (const Box& tbx, - Array4<Real const> const& tmp, - Array4<Real > const& dst, - int scomp, int dcomp, int ncomp) -{ - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); - Dim3 slen_local = slen; - AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, - { - Real d = 0.0; - - 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) - Real sss = sx[ix]*sy[iy]*sz[iz]; -#else - Real sss = sx[ix]*sz[iy]; -#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) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else - d += sss*( tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#endif - } - } - } - - dst(i,j,k,dcomp+n) = d; - }); -} - -#else - -void -BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - BL_PROFILE("BilinearFilter::ApplyStencil()"); - ncomp = std::min(ncomp, srcmf.nComp()); -#ifdef _OPENMP -#pragma omp parallel -#endif - { - FArrayBox tmpfab; - for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ - const auto& srcfab = srcmf[mfi]; - auto& dstfab = dstmf[mfi]; - const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - // tmpfab has enough ghost cells for the stencil - tmpfab.resize(gbx,ncomp); - tmpfab.setVal(0.0, gbx, 0, ncomp); - // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcfab.box(); - tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - // Apply filter - Filter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); - } - } -} - -void BilinearFilter::Filter (const Box& tbx, - Array4<Real const> const& tmp, - Array4<Real > const& dst, - int scomp, int dcomp, int ncomp) -{ - const auto lo = amrex::lbound(tbx); - const auto hi = amrex::ubound(tbx); - // tmp and dst are of type Array4 (Fortran ordering) - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); - for (int n = 0; n < ncomp; ++n) { - // Set dst value to 0. - for (int k = lo.z; k <= hi.z; ++k) { - for (int j = lo.y; j <= hi.y; ++j) { - for (int i = lo.x; i <= hi.x; ++i) { - dst(i,j,k,dcomp+n) = 0.0; - } - } - } - // 3 nested loop on 3D stencil - 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) - Real sss = sx[ix]*sy[iy]*sz[iz]; -#else - Real sss = sx[ix]*sz[iy]; -#endif - // 3 nested loop on 3D array - for (int k = lo.z; k <= hi.z; ++k) { - for (int j = lo.y; j <= hi.y; ++j) { - AMREX_PRAGMA_SIMD - for (int i = lo.x; i <= hi.x; ++i) { -#if (AMREX_SPACEDIM == 3) - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) - +tmp(i+ix,j-iy,k-iz,scomp+n) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#endif - } - } - } - } - } - } - } -} - -#endif diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H new file mode 100644 index 000000000..eaa0498c9 --- /dev/null +++ b/Source/Filter/Filter.H @@ -0,0 +1,43 @@ +#include <AMReX_MultiFab.H> +#include <AMReX_CudaContainers.H> + +#ifndef WARPX_FILTER_H_ +#define WARPX_FILTER_H_ + +class Filter +{ +public: + Filter () = default; + + // Apply stencil on MultiFab. + // Guard cells are handled inside this function + 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, + int scomp=0, int dcomp=0, int ncomp=10000); + + // public for cuda + void DoFilter(const amrex::Box& tbx, + amrex::Array4<amrex::Real const> const& tmp, + amrex::Array4<amrex::Real > const& dst, + int scomp, int dcomp, int ncomp); + + // In 2D, stencil_length_each_dir = {length(stencil_x), length(stencil_z)} + amrex::IntVect stencil_length_each_dir; + +protected: + // Stencil along each direction. + // in 2D, stencil_y is not initialized. + amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z; + // Length of each stencil. + // In 2D, slen = {length(stencil_x), length(stencil_z), 1} + amrex::Dim3 slen; + +private: + +}; +#endif // #ifndef WARPX_FILTER_H_ diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp new file mode 100644 index 000000000..f8a2dd6c2 --- /dev/null +++ b/Source/Filter/Filter.cpp @@ -0,0 +1,257 @@ +#include <WarpX.H> +#include <Filter.H> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +#ifdef AMREX_USE_CUDA + +/* \brief Apply stencil on MultiFab (GPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(MultiFab)"); + ncomp = std::min(ncomp, srcmf.nComp()); + + for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) + { + const auto& src = srcmf.array(mfi); + const auto& dst = dstmf.array(mfi); + const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcmf[mfi].box(); + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); + } +} + +/* \brief Apply stencil on FArrayBox (GPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + const auto& src = srcfab.array(); + const auto& dst = dstfab.array(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); +} + +/* \brief Apply stencil (2D/3D, CPU/GPU) + */ +void Filter::DoFilter (const Box& tbx, + Array4<Real const> const& tmp, + Array4<Real > const& dst, + int scomp, int dcomp, int ncomp) +{ + amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); + amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); + amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + Dim3 slen_local = slen; + AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, + { + Real d = 0.0; + + 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) + Real sss = sx[ix]*sy[iy]*sz[iz]; +#else + Real sss = sx[ix]*sz[iy]; +#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) + +tmp(i-ix,j+iy,k-iz,scomp+n) + +tmp(i+ix,j+iy,k-iz,scomp+n) + +tmp(i-ix,j-iy,k+iz,scomp+n) + +tmp(i+ix,j-iy,k+iz,scomp+n) + +tmp(i-ix,j+iy,k+iz,scomp+n) + +tmp(i+ix,j+iy,k+iz,scomp+n)); +#else + d += sss*( tmp(i-ix,j-iy,k,scomp+n) + +tmp(i+ix,j-iy,k,scomp+n) + +tmp(i-ix,j+iy,k,scomp+n) + +tmp(i+ix,j+iy,k,scomp+n)); +#endif + } + } + } + + dst(i,j,k,dcomp+n) = d; + }); +} + +#else + +/* \brief Apply stencil on MultiFab (CPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil()"); + ncomp = std::min(ncomp, srcmf.nComp()); +#ifdef _OPENMP +#pragma omp parallel +#endif + { + FArrayBox tmpfab; + for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ + const auto& srcfab = srcmf[mfi]; + auto& dstfab = dstmf[mfi]; + const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); + } + } +} + +/* \brief Apply stencil on FArrayBox (CPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + FArrayBox tmpfab; + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); +} + +void Filter::DoFilter (const Box& tbx, + Array4<Real const> const& tmp, + Array4<Real > const& dst, + int scomp, int dcomp, int ncomp) +{ + const auto lo = amrex::lbound(tbx); + const auto hi = amrex::ubound(tbx); + // tmp and dst are of type Array4 (Fortran ordering) + amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); + amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); + amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + for (int n = 0; n < ncomp; ++n) { + // Set dst value to 0. + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + dst(i,j,k,dcomp+n) = 0.0; + } + } + } + // 3 nested loop on 3D stencil + 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) + Real sss = sx[ix]*sy[iy]*sz[iz]; +#else + Real sss = sx[ix]*sz[iy]; +#endif + // 3 nested loop on 3D array + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + AMREX_PRAGMA_SIMD + for (int i = lo.x; i <= hi.x; ++i) { +#if (AMREX_SPACEDIM == 3) + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) + +tmp(i+ix,j-iy,k-iz,scomp+n) + +tmp(i-ix,j+iy,k-iz,scomp+n) + +tmp(i+ix,j+iy,k-iz,scomp+n) + +tmp(i-ix,j-iy,k+iz,scomp+n) + +tmp(i+ix,j-iy,k+iz,scomp+n) + +tmp(i-ix,j+iy,k+iz,scomp+n) + +tmp(i+ix,j+iy,k+iz,scomp+n)); +#else + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) + +tmp(i+ix,j-iy,k,scomp+n) + +tmp(i-ix,j+iy,k,scomp+n) + +tmp(i+ix,j+iy,k,scomp+n)); +#endif + } + } + } + } + } + } + } +} + +#endif // #ifdef AMREX_USE_CUDA diff --git a/Source/Filter/Make.package b/Source/Filter/Make.package index 36e74bfb0..8b8e9b50b 100644 --- a/Source/Filter/Make.package +++ b/Source/Filter/Make.package @@ -1,5 +1,9 @@ +CEXE_headers += Filter.H +CEXE_sources += Filter.cpp CEXE_sources += BilinearFilter.cpp CEXE_headers += BilinearFilter.H +CEXE_sources += NCIGodfreyFilter.cpp +CEXE_headers += NCIGodfreyFilter.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Filter VPATH_LOCATIONS += $(WARPX_HOME)/Source/Filter diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H new file mode 100644 index 000000000..a53039dfa --- /dev/null +++ b/Source/Filter/NCIGodfreyFilter.H @@ -0,0 +1,30 @@ +#include <Filter.H> + +#ifndef WARPX_GODFREY_FILTER_H_ +#define WARPX_GODFREY_FILTER_H_ + +enum class godfrey_coeff_set { Ex_Ey_Bz=0, Bx_By_Ez=1 }; + +class NCIGodfreyFilter : public Filter +{ +public: + + NCIGodfreyFilter () = default; + + NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_); + + void ComputeStencils(); + + void getGodfreyCoeffs(godfrey_coeff_set coeff_set_in); + + static constexpr int stencil_width = 4; + +private: + + godfrey_coeff_set coeff_set; + amrex::Real cdtodz; + int l_lower_order_in_v; + +}; + +#endif // #ifndef WARPX_GODFREY_FILTER_H_ diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp new file mode 100644 index 000000000..3f589260a --- /dev/null +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -0,0 +1,78 @@ +#include <WarpX.H> +#include <NCIGodfreyFilter.H> +#include <NCIGodfreyTables.H> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_){ + // Store parameters into class data members + coeff_set = coeff_set_; + cdtodz = cdtodz_; + l_lower_order_in_v = l_lower_order_in_v_; + // NCI Godfrey filter has fixed size, and is applied along z only. +#if (AMREX_SPACEDIM == 3) + stencil_length_each_dir = {1,1,5}; + slen = {1,1,5}; +#else + stencil_length_each_dir = {1,5}; + slen = {1,5,1}; +#endif +} + +void NCIGodfreyFilter::ComputeStencils(){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( +#if ( AMREX_SPACEDIM == 3 ) + slen.z==5, +#else + slen.y==5, +#endif + "ERROR: NCI filter requires 5 points in z"); + 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); + index = max(index, 0); + Real weight_right = cdtodz - index/tab_length; + 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] + + 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] + + weight_right*table_nci_godfrey_Bx_By_Ez[index+1][i]; + } + } + + // Compute stencil_z + stencil_z.resize( 5 ); + stencil_z[0] = (256 + 128*prestencil[0] + 96*prestencil[1] + 80*prestencil[2] + 70*prestencil[3]) / 256; + stencil_z[1] = -( 64*prestencil[0] + 64*prestencil[1] + 60*prestencil[2] + 56*prestencil[3]) / 256; + stencil_z[2] = ( 16*prestencil[1] + 24*prestencil[2] + 28*prestencil[3]) / 256; + 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, + // so only 1 coeff, equal to 1) + stencil_x.resize(1); + stencil_x[0] = 1.; +#if (AMREX_SPACEDIM == 3) + stencil_y.resize(1); + stencil_y[0] = 1.; +#endif + + // 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/FortranInterface/WarpX_f.F90 b/Source/FortranInterface/WarpX_f.F90 index ccd6dd48a..d7f7a8053 100644 --- a/Source/FortranInterface/WarpX_f.F90 +++ b/Source/FortranInterface/WarpX_f.F90 @@ -186,6 +186,41 @@ contains end subroutine warpx_compute_dive_2d + subroutine warpx_compute_dive_rz (lo, hi, dive, dlo, dhi, & + Ex, xlo, xhi, Ey, ylo, yhi, Ez, zlo, zhi, dx, rmin) & + bind(c, name='warpx_compute_dive_rz') + integer, intent(in) :: lo(2),hi(2),dlo(2),dhi(2),xlo(2),xhi(2),ylo(2),yhi(2),zlo(2),zhi(2) + real(amrex_real), intent(in) :: dx(3), rmin + real(amrex_real), intent(in ) :: Ex (xlo(1):xhi(1),xlo(2):xhi(2)) + real(amrex_real), intent(in ) :: Ey (ylo(1):yhi(1),ylo(2):yhi(2)) + real(amrex_real), intent(in ) :: Ez (zlo(1):zhi(1),zlo(2):zhi(2)) + real(amrex_real), intent(inout) :: dive(dlo(1):dhi(1),dlo(2):dhi(2)) + + integer :: i,k + real(amrex_real) :: dxinv(3) + real(amrex_real) :: ru, rd + + dxinv = 1.d0/dx + + do k = lo(2), hi(2) + do i = lo(1), hi(1) + if (i == 0 .and. rmin == 0.) then + ! the bulk equation diverges on axis + ! (due to the 1/r terms). The following expressions regularize + ! these divergences. + dive(i,k) = 4.*dxinv(1) * Ex(i,k) & + + dxinv(3) * (Ez(i,k) - Ez(i,k-1)) + else + ru = 1.d0 + 0.5d0/(rmin*dxinv(1) + i) + rd = 1.d0 - 0.5d0/(rmin*dxinv(1) + i) + dive(i,k) = dxinv(1) * (ru*Ex(i,k) - rd*Ex(i-1,k)) & + + dxinv(3) * (Ez(i,k) - Ez(i,k-1)) + end if + end do + end do + end subroutine warpx_compute_dive_rz + + subroutine warpx_sync_current_2d (lo, hi, crse, clo, chi, fine, flo, fhi, dir) & bind(c, name='warpx_sync_current_2d') integer, intent(in) :: lo(2), hi(2), flo(2), fhi(2), clo(2), chi(2), dir diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index f246f9f54..1053ace89 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -39,16 +39,9 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d -#define WRPX_LORENTZ_TRANSFORM_Z warpx_lorentz_transform_z -#define WRPX_COPY_SLICE warpx_copy_slice_3d -#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs -#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_3d - - #elif (AMREX_SPACEDIM == 2) #define WRPX_COMPUTE_DIVB warpx_compute_divb_2d -#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d #define WRPX_SYNC_CURRENT warpx_sync_current_2d #define WRPX_SYNC_RHO warpx_sync_rho_2d @@ -69,10 +62,11 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_2d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_2d -#define WRPX_LORENTZ_TRANSFORM_Z warpx_lorentz_transform_z -#define WRPX_COPY_SLICE warpx_copy_slice_2d -#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs -#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_2d +#ifdef WARPX_RZ +#define WRPX_COMPUTE_DIVE warpx_compute_dive_rz +#else +#define WRPX_COMPUTE_DIVE warpx_compute_dive_2d +#endif #endif @@ -87,11 +81,6 @@ extern "C" amrex_real* xpold, amrex_real* ypold, amrex_real* zpold, amrex_real* uxpold, amrex_real* uypold, amrex_real* uzpold); - void WRPX_COPY_SLICE(const int* lo, const int* hi, - const amrex_real* tmp, const int* tlo, const int* thi, - amrex_real* buf, const int* blo, const int* bhi, - const int* ncomp, const int* i_boost, const int* i_lab); - // Charge deposition void warpx_charge_deposition(amrex::Real* rho, const long* np, const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, const amrex::Real* w, @@ -102,6 +91,12 @@ extern "C" const long* nox, const long* noy,const long* noz, const long* lvect, const long* charge_depo_algo); + // Charge deposition finalize for RZ + void warpx_charge_deposition_rz_volume_scaling( + amrex::Real* rho, const long* rho_ng, const int* rho_ntot, + const amrex::Real* rmin, + const amrex::Real* dr); + // Current deposition void warpx_current_deposition( amrex::Real* jx, const long* jx_ng, const int* jx_ntot, @@ -165,12 +160,6 @@ extern "C" const amrex::Real* charge, const amrex::Real* mass, const amrex::Real* dt, const long* particle_pusher_algo); - // Particle pusher (position) - - void warpx_particle_pusher_positions(const long* np, - amrex::Real* xp, amrex::Real* yp, amrex::Real* zp, - amrex::Real* uxp, amrex::Real* uyp, amrex::Real* uzp, amrex::Real* gaminv, - const amrex::Real* dt); // Laser pusher @@ -344,10 +333,6 @@ extern "C" const amrex::Real* dt, const amrex::Real* prob_lo, const amrex::Real* prob_hi); - void WRPX_LORENTZ_TRANSFORM_Z(amrex::Real* data, const int* dlo, const int* dhi, - const int* tlo, const int* thi, - const amrex::Real* gamma_boost, const amrex::Real* beta_boost); - // These functions are used to evolve E and B in the PML void WRPX_COMPUTE_DIVB (const int* lo, const int* hi, @@ -362,7 +347,11 @@ extern "C" const BL_FORT_FAB_ARG_ANYD(ex), const BL_FORT_FAB_ARG_ANYD(ey), const BL_FORT_FAB_ARG_ANYD(ez), - const amrex::Real* dx); + const amrex::Real* dx +#ifdef WARPX_RZ + ,const amrex::Real* rmin +#endif + ); void WRPX_PUSH_PML_BVEC(const int* xlo, const int* xhi, const int* ylo, const int* yhi, @@ -459,14 +448,6 @@ extern "C" const BL_FORT_FAB_ARG_ANYD(fine), const int* ncomp); - void WRPX_PXR_NCI_CORR_INIT(amrex::Real*, amrex::Real*, const int, - const amrex::Real, const int); - - void WRPX_PXR_GODFREY_FILTER (const int* lo, const int* hi, - amrex_real* fou, const int* olo, const int* ohi, - const amrex_real* fin, const int* ilo, const int* ihi, - const amrex_real* stencil, const int* nsten); - #ifdef WARPX_USE_PSATD void warpx_fft_mpi_init (int fcomm); void warpx_fft_domain_decomp (int* warpx_local_nz, int* warpx_local_z0, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index 151342ff5..12d541b08 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -18,7 +18,8 @@ #define WRPX_PXR_GETEB_ENERGY_CONSERVING geteb2drz_energy_conserving_generic #define WRPX_PXR_CURRENT_DEPOSITION depose_jrjtjz_generic_rz -#define WRPX_PXR_RZ_VOLUME_SCALING apply_rz_volume_scaling +#define WRPX_PXR_RZ_VOLUME_SCALING_RHO apply_rz_volume_scaling_rho +#define WRPX_PXR_RZ_VOLUME_SCALING_J apply_rz_volume_scaling_j #define WRPX_PXR_PUSH_BVEC pxrpush_emrz_bvec #define WRPX_PXR_PUSH_EVEC pxrpush_emrz_evec @@ -116,7 +117,7 @@ contains pxr_ll4symtry = ll4symtry .eq. 1 pxr_l_lower_order_in_v = l_lower_order_in_v .eq. 1 pxr_l_nodal = l_nodal .eq. 1 - + exg_nguards = ixyzmin - exg_lo eyg_nguards = ixyzmin - eyg_lo ezg_nguards = ixyzmin - ezg_lo @@ -217,6 +218,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN CALL depose_rho_vecHVv2_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& nxguard,nyguard,nzguard,lvect) + + ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN + CALL depose_rho_vecHVv2_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& + nxguard,nyguard,nzguard,lvect) + ELSE CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& nxguard,nyguard,nzguard,nox,noy,noz, & @@ -227,9 +233,15 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! Dimension 2 #elif (AMREX_SPACEDIM==2) +#ifdef WARPX_RZ + logical(pxr_logical) :: l_2drz = .TRUE._c_long +#else + logical(pxr_logical) :: l_2drz = .FALSE._c_long +#endif + CALL pxr_depose_rho_n_2dxz(rho,np,xp,yp,zp,w,q,xmin,zmin,dx,dz,nx,nz,& nxguard,nzguard,nox,noz, & - .TRUE._c_long, .FALSE._c_long, .FALSE._c_long, 0_c_long) + .TRUE._c_long, .FALSE._c_long, l_2drz, 0_c_long) #endif @@ -238,6 +250,46 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! _________________________________________________________________ !> !> @brief + !> Applies the inverse volume scaling for RZ charge deposition + !> + !> @details + !> The scaling is done for both single mode (FDTD) and + !> multi mode (spectral) (todo) + ! + !> @param[inout] rho charge array + !> @param[in] rmin tile grid minimum radius + !> @param[in] dr radial space discretization steps + !> @param[in] nx,ny,nz number of cells + !> @param[in] nxguard,nyguard,nzguard number of guard cells + !> + subroutine warpx_charge_deposition_rz_volume_scaling(rho,rho_ng,rho_ntot,rmin,dr) & + bind(C, name="warpx_charge_deposition_rz_volume_scaling") + + integer, intent(in) :: rho_ntot(AMREX_SPACEDIM) + integer(c_long), intent(in) :: rho_ng + real(amrex_real), intent(IN OUT):: rho(*) + real(amrex_real), intent(IN) :: rmin, dr + +#ifdef WARPX_RZ + integer(c_long) :: type_rz_depose = 1 +#endif + + ! Compute the number of valid cells and guard cells + integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) + rho_nvalid = rho_ntot - 2*rho_ng + rho_nguards = rho_ng + +#ifdef WARPX_RZ + CALL WRPX_PXR_RZ_VOLUME_SCALING_RHO( & + rho,rho_nguards,rho_nvalid, & + rmin,dr,type_rz_depose) +#endif + + end subroutine warpx_charge_deposition_rz_volume_scaling + + ! _________________________________________________________________ + !> + !> @brief !> Main subroutine for the current deposition !> !> @details @@ -340,8 +392,9 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif ! Compute the number of valid cells and guard cells integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) @@ -353,7 +406,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n jz_nguards = jz_ng #ifdef WARPX_RZ - CALL WRPX_PXR_RZ_VOLUME_SCALING( & + CALL WRPX_PXR_RZ_VOLUME_SCALING_J( & jx,jx_nguards,jx_nvalid, & jy,jy_nguards,jy_nvalid, & jz,jz_nguards,jz_nvalid, & @@ -413,7 +466,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n END SELECT !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) #elif (AMREX_SPACEDIM == 2) CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) @@ -477,38 +530,6 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n ! _________________________________________________________________ !> !> @brief - !> Main subroutine for the particle pusher of positions - !> - !> @param[in] np number of super-particles - !> @param[in] xp,yp,zp particle position arrays - !> @param[in] uxp,uyp,uzp normalized momentum in each direction - !> @param[in] gaminv particle Lorentz factors - !> @param[in] dt time step - !> @param[in] particle_pusher_algo Particle pusher algorithm - subroutine warpx_particle_pusher_positions(np,xp,yp,zp,uxp,uyp,uzp, & - gaminv,dt) & - bind(C, name="warpx_particle_pusher_positions") - - INTEGER(c_long), INTENT(IN) :: np - REAL(amrex_real),INTENT(INOUT) :: gaminv(np) - REAL(amrex_real),INTENT(INOUT) :: xp(np),yp(np),zp(np) - REAL(amrex_real),INTENT(INOUT) :: uxp(np),uyp(np),uzp(np) - REAL(amrex_real),INTENT(IN) :: dt - - CALL pxr_set_gamma(np,uxp,uyp,uzp,gaminv) - - !!!! --- push particle species positions a time step -#if (AMREX_SPACEDIM == 3) - CALL pxr_pushxyz(np,xp,yp,zp,uxp,uyp,uzp,gaminv,dt) -#elif (AMREX_SPACEDIM == 2) - CALL pxr_pushxz(np,xp,zp,uxp,uzp,gaminv,dt) -#endif - - end subroutine warpx_particle_pusher_positions - - ! _________________________________________________________________ - !> - !> @brief !> Main subroutine for the Maxwell solver (E field) !> !> @param[in] xlo, xhi, ylo, yhi, zlo, zhi lower and higher bounds diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index cd5802a55..f998e217e 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -293,6 +293,7 @@ public: amrex::Real z_rms; amrex::Real q_tot; long npart; + int do_symmetrize = 0; bool radially_weighted = true; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 52b5d8b61..f9642d1b6 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -286,6 +286,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.get("z_rms", z_rms); pp.get("q_tot", q_tot); pp.get("npart", npart); + pp.query("do_symmetrize", do_symmetrize); gaussian_beam = true; parseMomentum(pp); } diff --git a/Source/Initialization/PlasmaProfiles.cpp b/Source/Initialization/PlasmaProfiles.cpp index e3382db06..d9d207f7e 100644 --- a/Source/Initialization/PlasmaProfiles.cpp +++ b/Source/Initialization/PlasmaProfiles.cpp @@ -17,21 +17,22 @@ Real PredefinedDensityProfile::getDensity(Real x, Real y, Real z) const { /// plateau between linear upramp and downramp, and parab transverse profile /// Real PredefinedDensityProfile::ParabolicChannel(Real x, Real y, Real z) const { - // params = [ramp_up plateau ramp_down rc n0] - Real ramp_up = params[0]; - Real plateau = params[1]; - Real ramp_down = params[2]; - Real rc = params[3]; - Real n0 = params[4]; + // params = [z_start ramp_up plateau ramp_down rc n0] + Real z_start = params[0]; + Real ramp_up = params[1]; + Real plateau = params[2]; + Real ramp_down = params[3]; + Real rc = params[4]; + Real n0 = params[5]; Real n; Real kp = PhysConst::q_e/PhysConst::c*sqrt( n0/(PhysConst::m_e*PhysConst::ep0) ); - if (z>=0 and z<ramp_up ) { - n = z/ramp_up; - } else if (z>=ramp_up and z<ramp_up+plateau ) { + if ((z-z_start)>=0 and (z-z_start)<ramp_up ) { + n = (z-z_start)/ramp_up; + } else if ((z-z_start)>=ramp_up and (z-z_start)<ramp_up+plateau ) { n = 1; - } else if (z>=ramp_up+plateau and z<ramp_up+plateau+ramp_down) { - n = 1-(z-ramp_up-plateau)/ramp_down; + } else if ((z-z_start)>=ramp_up+plateau and (z-z_start)<ramp_up+plateau+ramp_down) { + n = 1-((z-z_start)-ramp_up-plateau)/ramp_down; } else { n = 0; } diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 23637ec97..d583b2b0f 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -7,6 +7,7 @@ #include <WarpX.H> #include <WarpX_f.H> #include <BilinearFilter.H> +#include <NCIGodfreyFilter.H> #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> @@ -161,8 +162,6 @@ WarpX::InitNCICorrector () { if (WarpX::use_fdtd_nci_corr) { - mypc->fdtd_nci_stencilz_ex.resize(max_level+1); - mypc->fdtd_nci_stencilz_by.resize(max_level+1); for (int lev = 0; lev <= max_level; ++lev) { const Geometry& gm = Geom(lev); @@ -174,10 +173,15 @@ WarpX::InitNCICorrector () dz = dx[1]; } cdtodz = PhysConst::c * dt[lev] / dz; - WRPX_PXR_NCI_CORR_INIT( (mypc->fdtd_nci_stencilz_ex)[lev].data(), - (mypc->fdtd_nci_stencilz_by)[lev].data(), - mypc->nstencilz_fdtd_nci_corr, cdtodz, - WarpX::l_lower_order_in_v); + + // Initialize Godfrey filters + // Same filter for fields Ex, Ey and Bz + nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz, cdtodz, WarpX::l_lower_order_in_v) ); + // Same filter for fields Bx, By and Ez + nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez, cdtodz, WarpX::l_lower_order_in_v) ); + // Compute Godfrey filters stencils + nci_godfrey_filter_exeybz[lev]->ComputeStencils(); + nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); } } } @@ -318,6 +322,7 @@ WarpX::InitLevelData (int lev, Real time) void WarpX::InitLevelDataFFT (int lev, Real time) { + Efield_fp_fft[lev][0]->setVal(0.0); Efield_fp_fft[lev][1]->setVal(0.0); Efield_fp_fft[lev][2]->setVal(0.0); @@ -342,6 +347,7 @@ WarpX::InitLevelDataFFT (int lev, Real time) current_cp_fft[lev][2]->setVal(0.0); rho_cp_fft[lev]->setVal(0.0); } + } #endif diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 5dd94d0c7..516e368ab 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -84,11 +84,19 @@ private: std::string field_function; // laser particle domain - amrex::RealBox prob_domain; - + amrex::RealBox laser_injection_box; + // 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 + // outside of the simulation domain and enters it. + void ContinuousInjection(const amrex::RealBox& injection_box) override; + // Update position of the antenna + void UpdateContinuousInjectionPosition(amrex::Real dt) override; }; #endif diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 2b56c3cfd..2f964b6f3 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -25,8 +25,9 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, { charge = 1.0; mass = std::numeric_limits<Real>::max(); - - ParmParse pp(laser_name); + do_boosted_frame_diags = 0; + + ParmParse pp(laser_name); // Parse the type of laser profile and set the corresponding flag `profile` std::string laser_type_s; @@ -49,6 +50,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, pp.query("pusher_algo", pusher_algo); pp.get("wavelength", wavelength); pp.get("e_max", e_max); + pp.query("do_continuous_injection", do_continuous_injection); if ( profile == laser_t::Gaussian ) { // Parse the properties of the Gaussian profile @@ -76,14 +78,14 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, parser.define(field_function); parser.registerVariables({"X","Y","t"}); - ParmParse pp("my_constants"); + ParmParse ppc("my_constants"); std::set<std::string> symbols = parser.symbols(); symbols.erase("X"); symbols.erase("Y"); symbols.erase("t"); // after removing variables, we are left with constants for (auto it = symbols.begin(); it != symbols.end(); ) { Real v; - if (pp.query(it->c_str(), v)) { + if (ppc.query(it->c_str(), v)) { parser.setConstant(*it, v); it = symbols.erase(it); } else { @@ -148,18 +150,96 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, u_Y = {0., 1., 0.}; #endif - prob_domain = Geometry::ProbDomain(); + laser_injection_box= Geom(0).ProbDomain(); { Vector<Real> lo, hi; if (pp.queryarr("prob_lo", lo)) { - prob_domain.setLo(lo); + laser_injection_box.setLo(lo); } if (pp.queryarr("prob_hi", hi)) { - prob_domain.setHi(hi); + laser_injection_box.setHi(hi); + } + } + + if (do_continuous_injection){ + // 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); +#if (AMREX_SPACEDIM==2) + windir[2*dir] = 1.0; +#else + windir[dir] = 1.0; +#endif + 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[2]-1)*(WarpX::boost_direction[2]-1) < 1.e-12, + "do_continous_injection for laser particle only works if " + + "warpx.boost_direction = z. TODO: all directions."); } } } +/* \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 + // outdated full problem domain at t=0. + + // Convert updated_position to Real* to use RealBox::contains(). +#if (AMREX_SPACEDIM == 3) + const Real* p_pos = updated_position.dataPtr(); +#else + const Real p_pos[2] = {updated_position[0], updated_position[2]}; +#endif + if ( injection_box.contains(p_pos) ){ + // Update laser_injection_box with current value + laser_injection_box = injection_box; + // Inject laser particles. LaserParticleContainer::InitData + // is called only once, when the antenna enters the simulation + // domain. + InitData(); + } +} + +/* \brief update position of the antenna if running in boosted frame. + * \param dt time step (level 0). + * The up-to-date antenna position is stored in updated_position. + */ +void +LaserParticleContainer::UpdateContinuousInjectionPosition(Real dt) +{ + int dir = WarpX::moving_window_dir; + if (do_continuous_injection and (WarpX::gamma_boost > 1)){ + // In boosted-frame simulations, the antenna has moved since the last + // call to this function, and injection position needs to be updated +#if ( AMREX_SPACEDIM == 3 ) + updated_position[dir] -= WarpX::beta_boost * + WarpX::boost_direction[dir] * PhysConst::c * dt; +#elif ( AMREX_SPACEDIM == 2 ) + // In 2D, dir=0 corresponds to x and dir=1 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for both 2D and 3D simulations. + updated_position[2*dir] -= WarpX::beta_boost * + WarpX::boost_direction[2*dir] * PhysConst::c * dt; +#endif + } +} + void LaserParticleContainer::InitData () { @@ -175,6 +255,13 @@ LaserParticleContainer::InitData (int lev) ComputeSpacing(lev, S_X, S_Y); ComputeWeightMobility(S_X, S_Y); + // 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){ + position = updated_position; + } + auto Transform = [&](int i, int j) -> Vector<Real>{ #if (AMREX_SPACEDIM == 3) return { position[0] + (S_X*(i+0.5))*u_X[0] + (S_Y*(j+0.5))*u_Y[0], @@ -210,8 +297,8 @@ LaserParticleContainer::InitData (int lev) plane_hi[1] = std::max(plane_hi[1], j); }; - const Real* prob_lo = prob_domain.lo(); - const Real* prob_hi = prob_domain.hi(); + const Real* prob_lo = laser_injection_box.lo(); + const Real* prob_hi = laser_injection_box.hi(); #if (AMREX_SPACEDIM == 3) compute_min_max(prob_lo[0], prob_lo[1], prob_lo[2]); compute_min_max(prob_hi[0], prob_lo[1], prob_lo[2]); @@ -272,7 +359,7 @@ LaserParticleContainer::InitData (int lev) #else const Real x[2] = {pos[0], pos[2]}; #endif - if (prob_domain.contains(x)) + if (laser_injection_box.contains(x)) { for (int k = 0; k<2; ++k) { particle_x.push_back(pos[0]); @@ -342,8 +429,6 @@ LaserParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); - auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 95e2b4ec0..c1ef54c33 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -41,7 +41,7 @@ endif CFLAGS += -fPIC FFLAGS += -fPIC F90FLAGS += -fPIC - USERSuffix = .Lib + USERSuffix := .Lib DEFINES += -DWARPX_USE_PY endif @@ -107,7 +107,7 @@ ifeq ($(USE_PSATD),TRUE) INCLUDE_LOCATIONS += $(FFTW_HOME)/include LIBRARY_LOCATIONS += $(FFTW_HOME)/lib endif - USERSuffix += .PSATD + USERSuffix := $(USERSuffix).PSATD DEFINES += -DWARPX_USE_PSATD ifeq ($(USE_CUDA),TRUE) DEFINES += -DFFTW -DCUDA_FFT=1 # PICSAR uses it @@ -119,7 +119,7 @@ ifeq ($(USE_PSATD),TRUE) endif ifeq ($(USE_RZ),TRUE) - USERSuffix += .RZ + USERSuffix := $(USERSuffix).RZ DEFINES += -DWARPX_RZ endif @@ -154,17 +154,23 @@ else SHARED_OPTION = -shared endif -installwarpx: libwarpx$(DIM)d.so - mv libwarpx$(DIM)d.so Python/pywarpx - cd Python; python setup.py install --with-libwarpx $(DIM) $(PYINSTALLOPTIONS) +ifeq ($(USE_RZ),TRUE) + PYDIM = rz +else + PYDIM = $(DIM)d +endif + +installwarpx: libwarpx$(PYDIM).so + mv libwarpx$(PYDIM).so Python/pywarpx + cd Python; python setup.py install --with-libwarpx $(PYDIM) $(PYINSTALLOPTIONS) -libwarpx$(DIM)d.a: $(objForExecs) +libwarpx$(PYDIM).a: $(objForExecs) @echo Making static library $@ ... $(SILENT) $(AR) -crv $@ $^ $(SILENT) $(RM) AMReX_buildInfo.cpp @echo SUCCESS -libwarpx$(DIM)d.so: $(objForExecs) +libwarpx$(PYDIM).so: $(objForExecs) @echo Making dynamic library $@ ... $(SILENT) $(CXX) $(SHARED_OPTION) -fPIC -o $@ $^ $(LDFLAGS) $(libraries) $(SILENT) $(RM) AMReX_buildInfo.cpp diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 348b31c8b..00dcb85d0 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -79,7 +79,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, 1, ng); const Real* dx = Geom(lev-1).CellSize(); - const int ref_ratio = refRatio(lev-1)[0]; + const int refinement_ratio = refRatio(lev-1)[0]; #ifdef _OPENMP #pragma omp parallel #endif @@ -89,7 +89,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(refinement_ratio).refine(refinement_ratio); // so that ccbx is coarsenable const FArrayBox& cxfab = dBx[mfi]; const FArrayBox& cyfab = dBy[mfi]; @@ -106,18 +106,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &refinement_ratio,&use_limiter); #else amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[0]), BL_TO_FORTRAN_ANYD(bfab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &refinement_ratio,&use_limiter); amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio,&use_limiter); + &refinement_ratio,&use_limiter); #endif for (int idim = 0; idim < 3; ++idim) @@ -153,7 +153,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, 1, ng); MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, 1, ng); - const int ref_ratio = refRatio(lev-1)[0]; + const int refinement_ratio = refRatio(lev-1)[0]; #ifdef _OPEMP #pragma omp parallel #endif @@ -163,7 +163,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(refinement_ratio).refine(refinement_ratio); // so that ccbx is coarsenable const FArrayBox& cxfab = dEx[mfi]; const FArrayBox& cyfab = dEy[mfi]; @@ -180,18 +180,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &refinement_ratio,&use_limiter); #else amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[0]), BL_TO_FORTRAN_ANYD(efab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &refinement_ratio,&use_limiter); amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio); + &refinement_ratio); #endif for (int idim = 0; idim < 3; ++idim) @@ -246,7 +246,7 @@ void WarpX::FillBoundaryE (int lev, PatchType patch_type) { if (patch_type == PatchType::fine) - { + { if (do_pml && pml[lev]->ok()) { pml[lev]->ExchangeE(patch_type, @@ -364,7 +364,7 @@ WarpX::SyncCurrent () current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& refinement_ratio = refRatio(lev-1); std::array<const MultiFab*,3> fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -372,7 +372,7 @@ WarpX::SyncCurrent () std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, refinement_ratio[0]); } Vector<Array<std::unique_ptr<MultiFab>,3> > j_fp(finest_level+1); @@ -524,10 +524,10 @@ WarpX::SyncCurrent () void WarpX::SyncCurrent (const std::array<const amrex::MultiFab*,3>& fine, const std::array< amrex::MultiFab*,3>& crse, - int ref_ratio) + int refinement_ratio) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine[0]->nGrowVect() + 1) /ref_ratio; + BL_ASSERT(refinement_ratio == 2); + const IntVect& ng = (fine[0]->nGrowVect() + 1) /refinement_ratio; #ifdef _OPEMP #pragma omp parallel @@ -539,7 +539,7 @@ WarpX::SyncCurrent (const std::array<const amrex::MultiFab*,3>& fine, for (MFIter mfi(*crse[idim],true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,refinement_ratio),1); ffab.resize(fbx); fbx &= (*fine[idim])[mfi].box(); ffab.setVal(0.0); @@ -564,8 +564,8 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { rhoc[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rhof[lev], *rhoc[lev], ref_ratio[0]); + const IntVect& refinement_ratio = refRatio(lev-1); + SyncRho(*rhof[lev], *rhoc[lev], refinement_ratio[0]); } Vector<std::unique_ptr<MultiFab> > rho_f_g(finest_level+1); @@ -673,10 +673,10 @@ WarpX::SyncRho (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rhof, * averaging the values of the charge density of the fine patch (on the same level). */ void -WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) +WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int refinement_ratio) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine.nGrowVect()+1)/ref_ratio; + BL_ASSERT(refinement_ratio == 2); + const IntVect& ng = (fine.nGrowVect()+1)/refinement_ratio; const int nc = fine.nComp(); #ifdef _OPEMP @@ -687,7 +687,7 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) for (MFIter mfi(crse,true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,refinement_ratio),1); ffab.resize(fbx, nc); fbx &= fine[mfi].box(); ffab.setVal(0.0); @@ -710,7 +710,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& refinement_ratio = refRatio(lev-1); std::array<const MultiFab*,3> fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -718,7 +718,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, refinement_ratio[0]); } void @@ -824,8 +824,8 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev) { if (rho_fp[lev]) { rho_cp[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rho_fp[lev], *rho_cp[lev], ref_ratio[0]); + const IntVect& refinement_ratio = refRatio(lev-1); + SyncRho(*rho_fp[lev], *rho_cp[lev], refinement_ratio[0]); } } diff --git a/Source/Parser/WarpXParser.cpp b/Source/Parser/WarpXParser.cpp index 8c800215f..3237086f2 100644 --- a/Source/Parser/WarpXParser.cpp +++ b/Source/Parser/WarpXParser.cpp @@ -1,4 +1,5 @@ +#include <algorithm> #include "WarpXParser.H" WarpXParser::WarpXParser (std::string const& func_body) @@ -12,6 +13,8 @@ WarpXParser::define (std::string const& func_body) clear(); m_expression = func_body; + m_expression.erase(std::remove(m_expression.begin(),m_expression.end(),'\n'), + m_expression.end()); std::string f = m_expression + "\n"; #ifdef _OPENMP diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index acbfe3390..f33095738 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -10,5 +10,7 @@ CEXE_headers += WarpXParticleContainer.H CEXE_headers += RigidInjectedParticleContainer.H CEXE_headers += PhysicalParticleContainer.H +include $(WARPX_HOME)/Source/Particles/Pusher/Make.package + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index f3fd522a9..869126fef 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -129,9 +129,7 @@ public: void Checkpoint (const std::string& dir) const; - void WritePlotFile( const std::string& dir, - const amrex::Vector<int>& real_flags, - const amrex::Vector<std::string>& real_names) const; + void WritePlotFile (const std::string& dir) const; void Restart (const std::string& dir); @@ -156,6 +154,10 @@ public: int nSpecies() const {return nspecies;} + int nSpeciesBoostedFrameDiags() const {return nspecies_boosted_frame_diags;} + int mapSpeciesBoostedFrameDiags(int i) const {return map_species_boosted_frame_diags[i];} + int doBoostedFrameDiags() const {return do_boosted_frame_diags;} + int nSpeciesDepositOnMainGrid () const { int r = 0; for (int i : deposit_on_main_grid) { @@ -169,17 +171,14 @@ public: const amrex::Real z_old, const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, amrex::Vector<WarpXParticleContainer::DiagnosticParticleData>& parts) const; - - // - // Parameters for the Cherenkov corrector in the FDTD solver. - // Both stencils are calculated ar runtime. - // - // Number of coefficients for the stencil of the NCI corrector. - // The stencil is applied in the z direction only. - static constexpr int nstencilz_fdtd_nci_corr=5; - amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_ex; - amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_by; + // 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; + // Update injection position for continuously-injected species. + void UpdateContinuousInjectionPosition(amrex::Real dt) const; + int doContinuousInjection() const; std::vector<std::string> GetSpeciesNames() const { return species_names; } @@ -207,6 +206,13 @@ private: void ReadParameters (); + // Number of species dumped in BoostedFrameDiagnostics + int nspecies_boosted_frame_diags = 0; + // 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; + // runtime parameters int nlasers = 0; int nspecies = 1; // physical particles only. nspecies+nlasers == allcontainers.size(). diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index a4df1f83a..9d39ec2f9 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -8,8 +8,6 @@ using namespace amrex; -constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; - MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) { ReadParameters(); @@ -31,16 +29,31 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + // Compute the number of species for which lab-frame data is dumped + // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer + // particle IDs in map_species_lab_diags. + map_species_boosted_frame_diags.resize(nspecies); + nspecies_boosted_frame_diags = 0; + for (int i=0; i<nspecies; i++){ + auto& pc = allcontainers[i]; + if (pc->do_boosted_frame_diags){ + map_species_boosted_frame_diags[nspecies_boosted_frame_diags] = i; + do_boosted_frame_diags = 1; + nspecies_boosted_frame_diags += 1; + } + } + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - for (int i = 0; i < nspecies + nlasers; ++i) + for (int i = 0; i < nspecies_boosted_frame_diags; ++i) { - allcontainers[i]->AddRealComp("xold"); - allcontainers[i]->AddRealComp("yold"); - allcontainers[i]->AddRealComp("zold"); - allcontainers[i]->AddRealComp("uxold"); - allcontainers[i]->AddRealComp("uyold"); - allcontainers[i]->AddRealComp("uzold"); + int is = map_species_boosted_frame_diags[i]; + allcontainers[is]->AddRealComp("xold"); + allcontainers[is]->AddRealComp("yold"); + allcontainers[is]->AddRealComp("zold"); + allcontainers[is]->AddRealComp("uxold"); + allcontainers[is]->AddRealComp("uyold"); + allcontainers[is]->AddRealComp("uzold"); } pc_tmp->AddRealComp("xold"); pc_tmp->AddRealComp("yold"); @@ -376,13 +389,25 @@ MultiParticleContainer BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - for (int i = 0; i < nspecies; ++i){ - WarpXParticleContainer* pc = allcontainers[i].get(); + // Loop over particle species + for (int i = 0; i < nspecies_boosted_frame_diags; ++i){ + int isp = map_species_boosted_frame_diags[i]; + WarpXParticleContainer* pc = allcontainers[isp].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; 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 + // 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 + // 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()); @@ -414,3 +439,48 @@ 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 + * calls virtual function ContinuousInjection. + */ +void +MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const +{ + for (int i=0; i<nspecies+nlasers; i++){ + auto& pc = allcontainers[i]; + if (pc->do_continuous_injection){ + pc->ContinuousInjection(injection_box); + } + } +} + +/* \brief Update position of continuous injection parameters. + * \param dt: simulation time step (level 0) + * All classes inherited from WarpXParticleContainer do not have + * a position to update (PhysicalParticleContainer does not do anything). + */ +void +MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const +{ + for (int i=0; i<nspecies+nlasers; i++){ + auto& pc = allcontainers[i]; + if (pc->do_continuous_injection){ + pc->UpdateContinuousInjectionPosition(dt); + } + } +} + +int +MultiParticleContainer::doContinuousInjection() const +{ + int warpx_do_continuous_injection = 0; + for (int i=0; i<nspecies+nlasers; i++){ + auto& pc = allcontainers[i]; + if (pc->do_continuous_injection){ + warpx_do_continuous_injection = 1; + } + } + return warpx_do_continuous_injection; +} diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 362683879..bd8cfb47e 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -94,15 +94,18 @@ public: void AddGaussianBeam(amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, - amrex::Real q_tot, long npart); + amrex::Real q_tot, long npart, int do_symmetrize); + + void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, + 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 t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; - bool injected = false; - + virtual void ConvertUnits (ConvertDirection convert_dir) override; + protected: std::string species_name; @@ -122,6 +125,9 @@ protected: int GetRefineFac(const amrex::Real x, const amrex::Real y, const amrex::Real z); std::unique_ptr<amrex::IArrayBox> m_refined_injection_mask = nullptr; + // Inject particles during the whole simulation + void ContinuousInjection(const amrex::RealBox& injection_box) override; + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 17e6d98d9..1f517fccb 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -24,7 +24,7 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -81,6 +81,41 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_backward_propagation", do_backward_propagation); 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 + // for this species. + pp.query("do_boosted_frame_diags", do_boosted_frame_diags); + + pp.query("plot_species", plot_species); + int do_user_plot_vars; + 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 + // boosted frame + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ + plot_flags.resize(PIdx::nattribs + 6, 1); + } else { + plot_flags.resize(PIdx::nattribs, 1); + } + } else { + // Set plot_flag to 0 for all attribs + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ + plot_flags.resize(PIdx::nattribs + 6, 0); + } else { + plot_flags.resize(PIdx::nattribs, 0); + } + // 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), + "plot_vars argument not in ParticleStringNames"); + plot_flags[ParticleStringNames::to_index.at(var)] = 1; + } + } + } } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) @@ -146,7 +181,8 @@ 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); RealBox containing_bx = geom.ProbDomain(); @@ -156,13 +192,15 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution<double> disty(y_m, y_rms); std::normal_distribution<double> distz(z_m, z_rms); - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - if (ParallelDescriptor::IOProcessor()) { - std::array<Real, 3> u; - Real weight; - for (long i = 0; i < npart; ++i) { + std::array<Real, 3> u; + Real weight; + // 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; + } + for (long i = 0; i < npart; ++i) { #if ( AMREX_SPACEDIM == 3 | WARPX_RZ) weight = q_tot/npart/charge; Real x = distx(mt); @@ -174,29 +212,27 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real y = 0.; Real z = distz(mt); #endif - if (plasma_injector->insideBounds(x, y, z)) { - plasma_injector->getMomentum(u, x, y, z); - if (WarpX::gamma_boost > 1.) { - MapParticletoBoostedFrame(x, y, z, u); - } - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - attribs[PIdx::w ] = weight; - - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) - { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } - - AddOneParticle(0, 0, 0, x, y, z, attribs); + if (plasma_injector->insideBounds(x, y, z)) { + plasma_injector->getMomentum(u, x, y, z); + if (do_symmetrize){ + std::array<Real, 3> u_tmp; + Real x_tmp, y_tmp; + // Add four particles to the beam: + // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) + for (int ix=0; ix<2; ix++){ + for (int iy=0; iy<2; iy++){ + u_tmp = u; + x_tmp = x*std::pow(-1,ix); + u_tmp[0] *= std::pow(-1,ix); + y_tmp = y*std::pow(-1,iy); + u_tmp[1] *= std::pow(-1,iy); + CheckAndAddParticle(x_tmp, y_tmp, z, + u_tmp, weight/4); + } + } + } else { + CheckAndAddParticle(x, y, z, u, weight); + } } } } @@ -204,6 +240,39 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, } void +PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, + std::array<Real, 3> u, + Real weight) +{ + std::array<Real,PIdx::nattribs> attribs; + attribs.fill(0.0); + + // update attribs with input arguments + if (WarpX::gamma_boost > 1.) { + MapParticletoBoostedFrame(x, y, z, u); + } + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + attribs[PIdx::w ] = weight; + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + // need to create old values + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } + // add particle + AddOneParticle(0, 0, 0, x, y, z, attribs); +} + +void PhysicalParticleContainer::AddParticles (int lev) { BL_PROFILE("PhysicalParticleContainer::AddParticles()"); @@ -228,7 +297,8 @@ PhysicalParticleContainer::AddParticles (int lev) plasma_injector->y_rms, plasma_injector->z_rms, plasma_injector->q_tot, - plasma_injector->npart); + plasma_injector->npart, + plasma_injector->do_symmetrize); return; @@ -361,7 +431,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -468,7 +538,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -602,7 +672,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -710,7 +780,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -800,7 +870,6 @@ FieldGatherES (const amrex::Vector<std::array<std::unique_ptr<amrex::MultiFab>, const auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; const long np = pti.numParticles(); - auto& attribs = pti.GetAttribs(); auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; @@ -990,9 +1059,6 @@ PhysicalParticleContainer::FieldGather (int lev, { const std::array<Real,3>& dx = WarpX::CellSize(lev); - // WarpX assumes the same number of guard cells for Ex, Ey, Ez, Bx, By, Bz - long ng = Ex.nGrow(); - BL_ASSERT(OnSameGrids(lev,Ex)); MultiFab* cost = WarpX::getCosts(lev); @@ -1101,8 +1167,9 @@ PhysicalParticleContainer::Evolve (int lev, const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); - const auto& mypc = WarpX::GetInstance().GetPartContainer(); - const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + // 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; BL_ASSERT(OnSameGrids(lev,jx)); @@ -1134,7 +1201,7 @@ PhysicalParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); + const Box& box = pti.validbox(); auto& attribs = pti.GetAttribs(); @@ -1151,7 +1218,7 @@ PhysicalParticleContainer::Evolve (int lev, const long np = pti.numParticles(); - // Data on the grid + // Data on the grid FArrayBox const* exfab = &(Ex[pti]); FArrayBox const* eyfab = &(Ey[pti]); FArrayBox const* ezfab = &(Ez[pti]); @@ -1159,6 +1226,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* byfab = &(By[pti]); FArrayBox const* bzfab = &(Bz[pti]); + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1170,54 +1238,43 @@ PhysicalParticleContainer::Evolve (int lev, static_cast<int>(WarpX::noz)}); #endif - // both 2d and 3d + // Filter Ex (Both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD(Ex[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); + // Update exfab reference exfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD(Ez[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); ezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD(By[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); byfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD(Ey[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); eyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD(Bx[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); bxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD(Bz[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); bzfab = &filtered_Bz; #endif } @@ -1393,53 +1450,43 @@ PhysicalParticleContainer::Evolve (int lev, static_cast<int>(WarpX::noz)}); #endif - // both 2d and 3d + // Filter Ex (both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD((*cEx)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); + // Update exfab reference cexfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD((*cEz)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); cezfab = &filtered_Ez; + + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD((*cBy)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); cbyfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD((*cEy)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + 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)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD((*cBx)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + 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)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD((*cBz)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); cbzfab = &filtered_Bz; #endif } @@ -1683,7 +1730,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); @@ -1828,6 +1875,8 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real // Note the the slice should always move in the negative boost direction. AMREX_ALWAYS_ASSERT(z_new < z_old); + AMREX_ALWAYS_ASSERT(do_boosted_frame_diags == 1); + const int nlevs = std::max(0, finestLevel()+1); // we figure out a box for coarse-grained rejection. If the RealBox corresponding to a @@ -2004,3 +2053,14 @@ int PhysicalParticleContainer::GetRefineFac(const Real x, const Real y, const Re return ref_fac; } + +/* \brief Inject particles during the simulation + * \param injection_box: domain where particles should be injected. + */ +void +PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) +{ + // Inject plasma on level 0. Paticles will be redistributed. + const int lev=0; + AddPlasma(lev, injection_box); +} diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H new file mode 100644 index 000000000..42c61343e --- /dev/null +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -0,0 +1,76 @@ +#ifndef WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ +#define WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ + +#include <limits> +#include <WarpXParticleContainer.H> +#include <AMReX_REAL.H> + +#ifndef WARPX_RZ + +/* \brief Extract the particle's coordinates from the ParticleType struct `p`, + * and stores them in the variables `x`, `y`, `z`. */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void GetPosition( + amrex::Real& x, amrex::Real& y, amrex::Real& z, + const WarpXParticleContainer::ParticleType& p) +{ +#if (AMREX_SPACEDIM==3) + x = p.pos(0); + y = p.pos(1); + z = p.pos(2); +#else + x = p.pos(0); + y = std::numeric_limits<amrex::Real>::quiet_NaN(); + z = p.pos(1); +#endif +} + +/* \brief Set the particle's coordinates in the ParticleType struct `p`, + * from their values in the variables `x`, `y`, `z`. */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void SetPosition( + WarpXParticleContainer::ParticleType& p, + const amrex::Real x, const amrex::Real y, const amrex::Real z) +{ +#if (AMREX_SPACEDIM==3) + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; +#else + p.pos(0) = x; + p.pos(1) = z; +#endif +} + +# else // if WARPX_RZ is True + +/* \brief Extract the particle's coordinates from `theta` and the attributes + * of the ParticleType struct `p` (which contains the radius), + * and store them in the variables `x`, `y`, `z` */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void GetCartesianPositionFromCylindrical( + amrex::Real& x, amrex::Real& y, amrex::Real& z, + const WarpXParticleContainer::ParticleType& p, const amrex::Real theta) +{ + const amrex::Real r = p.pos(0); + x = r*std::cos(theta); + y = r*std::sin(theta); + z = p.pos(1); +} + +/* \brief Set the particle's cylindrical coordinates by setting `theta` + * and the attributes of the ParticleType struct `p` (which stores the radius), + * from the values of `x`, `y`, `z` */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void SetCylindricalPositionFromCartesian( + WarpXParticleContainer::ParticleType& p, amrex::Real& theta, + const amrex::Real x, const amrex::Real y, const amrex::Real z ) +{ + theta = std::atan2(y, x); + p.pos(0) = std::sqrt(x*x + y*y); + p.pos(1) = z; +} + +#endif // WARPX_RZ + +#endif // WARPX_PARTICLES_PUSHER_GETANDSETPOSITION_H_ diff --git a/Source/Particles/Pusher/Make.package b/Source/Particles/Pusher/Make.package new file mode 100644 index 000000000..8c8e77905 --- /dev/null +++ b/Source/Particles/Pusher/Make.package @@ -0,0 +1,4 @@ +CEXE_headers += GetAndSetPosition.H +CEXE_headers += UpdatePosition.H +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Pusher diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H new file mode 100644 index 000000000..0a4f579f4 --- /dev/null +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -0,0 +1,30 @@ +#ifndef WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ +#define WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ + +#include <AMReX_FArrayBox.H> +#include <WarpXConst.H> +#include <AMReX_REAL.H> + +/* \brief Push the particle's positions over one timestep, + * given the value of its momenta `ux`, `uy`, `uz` */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdatePosition( + amrex::Real& x, amrex::Real& y, amrex::Real& z, + const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, + const amrex::Real dt ) +{ + + constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); + + // Compute inverse Lorentz factor + const amrex::Real inv_gamma = 1./std::sqrt(1. + (ux*ux + uy*uy + uz*uz)*inv_c2); + // Update positions over one time step + x += ux * inv_gamma * dt; +#if (AMREX_SPACEDIM == 3) || (defined WARPX_RZ) // RZ pushes particles in 3D + y += uy * inv_gamma * dt; +#endif + z += uz * inv_gamma * dt; + +} + +#endif // WARPX_PARTICLES_PUSHER_UPDATEPOSITION_H_ diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index d3a69f5b0..0b27a2f2f 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -57,6 +57,10 @@ public: const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; + virtual void ReadHeader (std::istream& is) override; + + virtual void WriteHeader (std::ostream& os) const override; + private: // User input quantities diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a5acca281..2a3e8dd0d 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -77,7 +77,6 @@ RigidInjectedParticleContainer::RemapParticles() for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - auto& attribs = pti.GetAttribs(); auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; @@ -225,7 +224,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 275554cd8..1e3dfb2ae 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -6,6 +6,8 @@ #include <AMReX_Particles.H> #include <AMReX_AmrCore.H> +enum struct ConvertDirection{WarpX_to_SI, SI_to_WarpX}; + struct PIdx { enum { // Particle Attributes stored in amrex::ParticleContainer's struct of array @@ -59,7 +61,6 @@ 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(); } @@ -85,7 +86,13 @@ class WarpXParticleContainer public: friend MultiParticleContainer; + // 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 + // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector<std::map<std::pair<int, int>, DiagnosticParticleData> >; WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); @@ -183,7 +190,18 @@ public: int thread_num, int 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. + // Current status: + // PhysicalParticleContainer: implemented. + // LaserParticleContainer: implemented. + // RigidInjectedParticleContainer: not implemented. + virtual void ContinuousInjection(const amrex::RealBox& injection_box) {} + // 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. @@ -207,9 +225,11 @@ public: amrex::Real x, amrex::Real y, amrex::Real z, std::array<amrex::Real,PIdx::nattribs>& attribs); - void ReadHeader (std::istream& is); + virtual void ReadHeader (std::istream& is); + + virtual void WriteHeader (std::ostream& os) const; - void WriteHeader (std::ostream& os) const; + virtual void ConvertUnits (ConvertDirection convert_dir){}; static void ReadParameters (); @@ -235,6 +255,8 @@ public: AddIntComp(comm); } + int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + protected: std::map<std::string, int> particle_comps; @@ -248,6 +270,14 @@ protected: static int do_not_push; + // 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 + // support all features allowed by direct injection. + int do_continuous_injection = 0; + + int do_boosted_frame_diags = 1; + amrex::Vector<amrex::FArrayBox> local_rho; amrex::Vector<amrex::FArrayBox> local_jx; amrex::Vector<amrex::FArrayBox> local_jy; @@ -255,9 +285,19 @@ protected: amrex::Vector<amrex::Cuda::ManagedDeviceVector<amrex::Real> > m_xp, m_yp, m_zp, m_giv; + // 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 + // to dump to file. + amrex::Vector<int> plot_flags; + // list of names of attributes to dump. + amrex::Vector<std::string> plot_vars; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; + }; #endif diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 2edd3c636..47d57294d 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -7,6 +7,10 @@ #include <WarpX_f.H> #include <WarpX.H> +// Import low-level single-particle kernels +#include <GetAndSetPosition.H> +#include <UpdatePosition.H> + using namespace amrex; int WarpXParticleContainer::do_not_push = 0; @@ -78,7 +82,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -86,9 +90,9 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["uxold"] = PIdx::nattribs+3; particle_comps["uyold"] = PIdx::nattribs+4; particle_comps["uzold"] = PIdx::nattribs+5; - + } - + // Initialize temporary local arrays for charge/current deposition int num_threads = 1; #ifdef _OPENMP @@ -174,7 +178,7 @@ WarpXParticleContainer::AddNParticles (int lev, int n, const Real* x, const Real* y, const Real* z, const Real* vx, const Real* vy, const Real* vz, int nattr, const Real* attr, int uniqueparticles, int id) -{ +{ BL_ASSERT(nattr == 1); const Real* weight = attr; @@ -230,15 +234,15 @@ WarpXParticleContainer::AddNParticles (int lev, #endif p.pos(1) = z[i]; #endif - - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x[i]); - particle_tile.push_back_real(particle_comps["yold"], y[i]); - particle_tile.push_back_real(particle_comps["zold"], z[i]); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + ptile.push_back_real(particle_comps["xold"], x[i]); + ptile.push_back_real(particle_comps["yold"], y[i]); + ptile.push_back_real(particle_comps["zold"], z[i]); } - + particle_tile.push_back(p); } @@ -249,14 +253,14 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); } - + for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) { #ifdef WARPX_RZ @@ -327,7 +331,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, jx_ptr = local_jx[thread_num].dataPtr(); jy_ptr = local_jy[thread_num].dataPtr(); jz_ptr = local_jz[thread_num].dataPtr(); - + local_jx[thread_num].setVal(0.0); local_jy[thread_num].setVal(0.0); local_jz[thread_num].setVal(0.0); @@ -426,7 +430,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, BL_PROFILE_VAR_STOP(blp_pxr_cd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); jx[pti].atomicAdd(local_jx[thread_num], tbx, tbx, 0, 0, 1); @@ -477,7 +481,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, auto jyntot = local_jy[thread_num].length(); auto jzntot = local_jz[thread_num].length(); #endif - + long ncrse = np - np_current; BL_PROFILE_VAR_START(blp_pxr_cd); if (j_is_nodal) { @@ -608,7 +612,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array<Real, 3>& xyzmin = xyzmin_tile; #ifdef AMREX_USE_GPU - data_ptr = (*rhomf)[pti].dataPtr(); + data_ptr = (*rhomf)[pti].dataPtr(icomp); auto rholen = (*rhomf)[pti].length(); #else tile_box.grow(ngRho); @@ -641,9 +645,14 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); (*rhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); @@ -660,7 +669,7 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, const std::array<Real,3>& cxyzmin_tile = WarpX::LowerCorner(ctilebox, lev-1); #ifdef AMREX_USE_GPU - data_ptr = (*crhomf)[pti].dataPtr(); + data_ptr = (*crhomf)[pti].dataPtr(icomp); auto rholen = (*crhomf)[pti].length(); #else tile_box = amrex::convert(ctilebox, IntVect::TheUnitVector()); @@ -696,9 +705,14 @@ WarpXParticleContainer::DepositCharge ( WarpXParIter& pti, RealVector& wp, &ngRho, &ngRho, &ngRho, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &cxyzmin_tile[0], &cdx[0]); +#endif BL_PROFILE_VAR_STOP(blp_pxr_chd); -#ifndef AMREX_USE_GPU +#ifndef AMREX_USE_GPU BL_PROFILE_VAR_START(blp_accumulate); (*crhomf)[pti].atomicAdd(local_rho[thread_num], tile_box, tile_box, 0, icomp, 1); @@ -723,7 +737,6 @@ WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, const auto& gm = m_gdb->Geom(lev); const auto& ba = m_gdb->ParticleBoxArray(lev); - const auto& dm = m_gdb->DistributionMap(lev); const Real* dx = gm.CellSize(); const Real* plo = gm.ProbLo(); @@ -793,36 +806,36 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #ifdef _OPENMP #pragma omp parallel -#endif { +#endif Cuda::ManagedDeviceVector<Real> xp, yp, zp; - FArrayBox local_rho; +#ifdef _OPENMP + FArrayBox rho_loc; +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = pti.validbox(); - auto& wp = pti.GetAttribs(PIdx::w); const long np = pti.numParticles(); pti.GetPosition(xp, yp, zp); - const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); - // Data on the grid Real* data_ptr; FArrayBox& rhofab = (*rho)[pti]; #ifdef _OPENMP + const std::array<Real,3>& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); const std::array<Real, 3>& xyzmin = xyzmin_tile; tile_box.grow(ng); - local_rho.resize(tile_box); - local_rho = 0.0; - data_ptr = local_rho.dataPtr(); - auto rholen = local_rho.length(); + rho_loc.resize(tile_box); + rho_loc = 0.0; + data_ptr = rho_loc.dataPtr(); + auto rholen = rho_loc.length(); #else + const Box& box = pti.validbox(); + const std::array<Real,3>& xyzmin_grid = WarpX::LowerCorner(box, lev); const std::array<Real, 3>& xyzmin = xyzmin_grid; data_ptr = rhofab.dataPtr(); auto rholen = rhofab.length(); @@ -852,12 +865,17 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) &dx[0], &dx[1], &dx[2], &nx, &ny, &nz, &nxg, &nyg, &nzg, &WarpX::nox,&WarpX::noy,&WarpX::noz, &lvect, &WarpX::charge_deposition_algo); +#ifdef WARPX_RZ + long ngRho = WarpX::nox; + warpx_charge_deposition_rz_volume_scaling( + data_ptr, &ngRho, rholen.getVect(), + &xyzmin[0], &dx[0]); +#endif #ifdef _OPENMP - rhofab.atomicAdd(local_rho); -#endif + rhofab.atomicAdd(rho_loc); } - +#endif } if (!local) rho->SumBoundary(gm.periodicity()); @@ -1007,8 +1025,6 @@ void WarpXParticleContainer::PushX (int lev, Real dt) { BL_PROFILE("WPC::PushX()"); - BL_PROFILE_VAR_NS("WPC::PushX::Copy", blp_copy); - BL_PROFILE_VAR_NS("WPC:PushX::Push", blp_pxr_pp); if (do_not_push) return; @@ -1018,47 +1034,42 @@ WarpXParticleContainer::PushX (int lev, Real dt) #pragma omp parallel #endif { - Cuda::ManagedDeviceVector<Real> xp, yp, zp, giv; for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { Real wt = amrex::second(); - auto& attribs = pti.GetAttribs(); - - auto& uxp = attribs[PIdx::ux]; - auto& uyp = attribs[PIdx::uy]; - auto& uzp = attribs[PIdx::uz]; - - const long np = pti.numParticles(); - - giv.resize(np); - - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); - pti.GetPosition(xp, yp, zp); - BL_PROFILE_VAR_STOP(blp_copy); - // // Particle Push // - BL_PROFILE_VAR_START(blp_pxr_pp); - warpx_particle_pusher_positions(&np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - giv.dataPtr(), &dt); - BL_PROFILE_VAR_STOP(blp_pxr_pp); - - // - // copy particle data back - // - BL_PROFILE_VAR_START(blp_copy); - pti.SetPosition(xp, yp, zp); - BL_PROFILE_VAR_STOP(blp_copy); + // Extract pointers to particle position and momenta, for this particle tile + // - positions are stored as an array of struct, in `ParticleType` + ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); + // - momenta are stored as a struct of array, in `attribs` + auto& attribs = pti.GetAttribs(); + Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); +#ifdef WARPX_RZ + Real* AMREX_RESTRICT theta = attribs[PIdx::theta].dataPtr(); +#endif + // Loop over the particles and update their position + amrex::ParallelFor( pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + ParticleType& p = pstructs[i]; // Particle object that gets updated + Real x, y, z; // Temporary variables +#ifndef WARPX_RZ + GetPosition( x, y, z, p ); // Initialize x, y, z + UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); + SetPosition( p, x, y, z ); // Update the object p +#else + // For WARPX_RZ, the particles are still pushed in 3D Cartesian + GetCartesianPositionFromCylindrical( x, y, z, p, theta[i] ); + UpdatePosition( x, y, z, ux[i], uy[i], uz[i], dt); + SetCylindricalPositionFromCartesian( p, theta[i], x, y, z ); +#endif + } + ); if (cost) { const Box& tbx = pti.tilebox(); @@ -1074,15 +1085,15 @@ WarpXParticleContainer::PushX (int lev, Real dt) } // This function is called in Redistribute, just after locate -void -WarpXParticleContainer::particlePostLocate(ParticleType& p, +void +WarpXParticleContainer::particlePostLocate(ParticleType& p, const ParticleLocData& pld, const int lev) { // Tag particle if goes to higher level. // It will be split later in the loop - if (pld.m_lev == lev+1 - and p.m_idata.id != NoSplitParticleID + if (pld.m_lev == lev+1 + and p.m_idata.id != NoSplitParticleID and p.m_idata.id >= 0) { p.m_idata.id = DoSplitParticleID; diff --git a/Source/Python/Make.package b/Source/Python/Make.package index 71bd4ebe8..746257abf 100644 --- a/Source/Python/Make.package +++ b/Source/Python/Make.package @@ -1,7 +1,7 @@ -ifeq ($(USE_PYTHON_MAIN),TRUE) - CEXE_sources += WarpXWrappers.cpp - CEXE_headers += WarpXWrappers.h -endif +#ifeq ($(USE_PYTHON_MAIN),TRUE) +# CEXE_sources += WarpXWrappers.cpp +# CEXE_headers += WarpXWrappers.h +#endif CEXE_sources += WarpXWrappers.cpp CEXE_sources += WarpX_py.cpp CEXE_headers += WarpXWrappers.h diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 16d7cd841..3c1a930b3 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -89,7 +89,7 @@ extern "C" void amrex_finalize (int finalize_mpi) { - amrex::Finalize(finalize_mpi); + amrex::Finalize(); } void warpx_init () diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index d8e2d2dab..cd335dbcf 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -4,6 +4,9 @@ CEXE_sources += WarpXTagging.cpp CEXE_sources += WarpXUtil.cpp CEXE_headers += WarpXConst.H CEXE_headers += WarpXUtil.H +CEXE_headers += WarpXAlgorithmSelection.H +CEXE_sources += WarpXAlgorithmSelection.cpp +CEXE_headers += NCIGodfreyTables.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H new file mode 100644 index 000000000..8cb105aa0 --- /dev/null +++ b/Source/Utils/NCIGodfreyTables.H @@ -0,0 +1,224 @@ +#include <AMReX_AmrCore.H> + +#ifndef WARPX_GODFREY_COEFF_TABLE_H_ +#define WARPX_GODFREY_COEFF_TABLE_H_ + +// Table width. This is related to the stencil length +const int tab_width = 4; +// table length. Each line correspond to 1 value of cdtodz +// (here 101 values). +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.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.65437,2.44244,-0.878945,0.0914107 + }; + +#endif // #ifndef WARPX_GODFREY_COEFF_TABLE_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H new file mode 100644 index 000000000..3fb23698a --- /dev/null +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -0,0 +1,57 @@ +#ifndef UTILS_WARPXALGORITHMSELECTION_H_ +#define UTILS_WARPXALGORITHMSELECTION_H_ + +#include <AMReX_ParmParse.H> +#include <string> + +struct MaxwellSolverAlgo { + // These numbers corresponds to the algorithm code in WarpX's + // `warpx_push_bvec` and `warpx_push_evec_f` + enum { + Yee = 0, + CKC = 1 + }; +}; + +struct ParticlePusherAlgo { + // These numbers correspond to the algorithm code in WarpX's + // `warpx_particle_pusher` + enum { + Boris = 0, + Vay = 1 + }; +}; + +struct CurrentDepositionAlgo { + enum { + // These numbers corresponds to the algorithm code in PICSAR's + // `depose_jxjyjz_generic` and `depose_jxjyjz_generic_2d` + Direct = 3, + DirectVectorized = 2, + EsirkepovNonOptimized = 1, + Esirkepov = 0 + }; +}; + +struct ChargeDepositionAlgo { + // These numbers corresponds to the algorithm code in WarpX's + // `warpx_charge_deposition` function + enum { + Vectorized = 0, + Standard = 1 + }; +}; + +struct GatheringAlgo { + // These numbers corresponds to the algorithm code in PICSAR's + // `geteb3d_energy_conserving_generic` function + enum { + Vectorized = 0, + Standard = 1 + }; +}; + +int +GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ); + +#endif // UTILS_WARPXALGORITHMSELECTION_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp new file mode 100644 index 000000000..2c8038ccd --- /dev/null +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -0,0 +1,95 @@ +#include <WarpXAlgorithmSelection.H> +#include <map> +#include <algorithm> +#include <cstring> + +// Define dictionary with correspondance between user-input strings, +// and corresponding integer for use inside the code (e.g. in PICSAR). + +const std::map<std::string, int> maxwell_solver_algo_to_int = { + {"yee", MaxwellSolverAlgo::Yee }, +#ifndef WARPX_RZ // Not available in RZ + {"ckc", MaxwellSolverAlgo::CKC }, +#endif + {"default", MaxwellSolverAlgo::Yee } +}; + +const std::map<std::string, int> particle_pusher_algo_to_int = { + {"boris", ParticlePusherAlgo::Boris }, + {"vay", ParticlePusherAlgo::Vay }, + {"default", ParticlePusherAlgo::Boris } +}; + +const std::map<std::string, int> current_deposition_algo_to_int = { + {"esirkepov", CurrentDepositionAlgo::Esirkepov }, + {"direct", CurrentDepositionAlgo::Direct }, +#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D + {"direct-vectorized", CurrentDepositionAlgo::DirectVectorized }, +#endif + {"default", CurrentDepositionAlgo::Esirkepov } +}; + +const std::map<std::string, int> charge_deposition_algo_to_int = { + {"standard", ChargeDepositionAlgo::Standard }, +#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D + {"vectorized", ChargeDepositionAlgo::Vectorized }, + {"default", ChargeDepositionAlgo::Vectorized } +#else + {"default", ChargeDepositionAlgo::Standard } +#endif +}; + +const std::map<std::string, int> gathering_algo_to_int = { + {"standard", GatheringAlgo::Standard }, +#ifndef AMREX_USE_GPU // Only available on CPU + {"vectorized", GatheringAlgo::Vectorized }, + {"default", GatheringAlgo::Vectorized } +#else + {"default", GatheringAlgo::Standard } +#endif +}; + + +int +GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ + + // Read user input ; use "default" if it is not found + std::string algo = "default"; + pp.query( pp_search_key, algo ); + // Convert to lower case + std::transform(algo.begin(), algo.end(), algo.begin(), ::tolower); + + // Pick the right dictionary + std::map<std::string, int> algo_to_int; + if (0 == std::strcmp(pp_search_key, "maxwell_fdtd_solver")) { + algo_to_int = maxwell_solver_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "particle_pusher")) { + algo_to_int = particle_pusher_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "current_deposition")) { + algo_to_int = current_deposition_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "charge_deposition")) { + algo_to_int = charge_deposition_algo_to_int; + } else if (0 == std::strcmp(pp_search_key, "field_gathering")) { + algo_to_int = gathering_algo_to_int; + } else { + std::string pp_search_string = pp_search_key; + amrex::Abort("Unknown algorithm type: " + pp_search_string); + } + + // Check if the user-input is a valid key for the dictionary + if (algo_to_int.count(algo) == 0){ + // Not a valid key ; print error message + std::string pp_search_string = pp_search_key; + std::string error_message = "Invalid string for algo." + pp_search_string + + ": " + algo + ".\nThe valid values are:\n"; + for ( const auto &valid_pair : algo_to_int ) { + if (valid_pair.first != "default"){ + error_message += " - " + valid_pair.first + "\n"; + } + } + amrex::Abort(error_message); + } + + // If the input is a valid key, return the value + return algo_to_int[algo]; +} diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 05e171f22..ae781f9aa 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -5,21 +5,21 @@ using namespace amrex; void -WarpX::UpdatePlasmaInjectionPosition (Real dt) +WarpX::UpdatePlasmaInjectionPosition (Real a_dt) { int dir = moving_window_dir; // Continuously inject plasma in new cells (by default only on level 0) - if (WarpX::do_plasma_injection and (WarpX::gamma_boost > 1)){ + if (WarpX::warpx_do_continuous_injection and (WarpX::gamma_boost > 1)){ // In boosted-frame simulations, the plasma has moved since the last // call to this function, and injection position needs to be updated current_injection_position -= WarpX::beta_boost * #if ( AMREX_SPACEDIM == 3 ) - WarpX::boost_direction[dir] * PhysConst::c * dt; + WarpX::boost_direction[dir] * PhysConst::c * a_dt; #elif ( AMREX_SPACEDIM == 2 ) // In 2D, dir=0 corresponds to x and dir=1 corresponds to z // This needs to be converted in order to index `boost_direction` // which has 3 components, for both 2D and 3D simulations. - WarpX::boost_direction[2*dir] * PhysConst::c * dt; + WarpX::boost_direction[2*dir] * PhysConst::c * a_dt; #endif } } @@ -33,7 +33,16 @@ WarpX::MoveWindow (bool move_j) // and of the plasma injection moving_window_x += moving_window_v * dt[0]; int dir = moving_window_dir; + + // Update warpx.current_injection_position + // PhysicalParticleContainer uses this injection position UpdatePlasmaInjectionPosition( dt[0] ); + if (WarpX::warpx_do_continuous_injection){ + // Update injection position for WarpXParticleContainer in mypc. + // Nothing to do for PhysicalParticleContainers + // For LaserParticleContainer, need to update the antenna position. + mypc->UpdateContinuousInjectionPosition( dt[0] ); + } // compute the number of cells to shift on the base level Real new_lo[AMREX_SPACEDIM]; @@ -53,8 +62,28 @@ WarpX::MoveWindow (bool move_j) } new_lo[dir] = current_lo[dir] + num_shift_base * cdx[dir]; new_hi[dir] = current_hi[dir] + num_shift_base * cdx[dir]; - RealBox new_box(new_lo, new_hi); - Geometry::ProbDomain(new_box); + + ResetProbDomain(RealBox(new_lo, new_hi)); + + // Moving slice coordinates - lo and hi - with moving window // + // slice box is modified only if slice diagnostics is initialized in input // + if ( slice_plot_int > 0 ) + { + Real new_slice_lo[AMREX_SPACEDIM]; + Real new_slice_hi[AMREX_SPACEDIM]; + const Real* current_slice_lo = slice_realbox.lo(); + const Real* current_slice_hi = slice_realbox.hi(); + for ( int i = 0; i < AMREX_SPACEDIM; i++) { + 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 - + 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]; + slice_realbox.setLo(new_slice_lo); + slice_realbox.setHi(new_slice_hi); + } int num_shift = num_shift_base; int num_shift_crse = num_shift; @@ -134,7 +163,7 @@ WarpX::MoveWindow (bool move_j) } // Continuously inject plasma in new cells (by default only on level 0) - if (WarpX::do_plasma_injection) { + if (WarpX::warpx_do_continuous_injection) { const int lev = 0; @@ -162,15 +191,11 @@ WarpX::MoveWindow (bool move_j) particleBox.setLo( dir, new_injection_position ); particleBox.setHi( dir, current_injection_position ); } - // Perform the injection of new particles in particleBox + if (particleBox.ok() and (current_injection_position != new_injection_position)){ - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast<PhysicalParticleContainer&>(pc); - ppc.AddPlasma(lev, particleBox); - } - // Update the injection position + // Performs continuous injection of all WarpXParticleContainer + // in mypc. + mypc->ContinuousInjection(particleBox); current_injection_position = new_injection_position; } } @@ -249,3 +274,14 @@ WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir) }); } } + +void +WarpX::ResetProbDomain (const RealBox& rb) +{ + Geometry::ResetDefaultProbDomain(rb); + for (int lev = 0; lev <= max_level; ++lev) { + Geometry g = Geom(lev); + g.ProbDomain(rb); + SetGeometry(lev, g); + } +} diff --git a/Source/Utils/WarpXTagging.cpp b/Source/Utils/WarpXTagging.cpp index 21acbefdc..8ea3211a3 100644 --- a/Source/Utils/WarpXTagging.cpp +++ b/Source/Utils/WarpXTagging.cpp @@ -9,7 +9,7 @@ using namespace amrex; void WarpX::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/) { - const Real* problo = Geometry::ProbLo(); + const Real* problo = Geom(lev).ProbLo(); const Real* dx = Geom(lev).CellSize(); #ifdef _OPENMP diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4ec7ebb51..a5ea6d75a 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -118,8 +118,7 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) const int lo_ind = bx.loVect()[1]; #endif // Check if box intersect with [zmin, zmax] - if ( (zmin>zmin_box && zmin<=zmax_box) || - (zmax>zmin_box && zmax<=zmax_box) ){ + if ( (zmax>zmin_box && zmin<=zmax_box) ){ Array4<Real> arr = mf[mfi].array(); // Set field to 0 between zmin and zmax ParallelFor(bx, diff --git a/Source/WarpX.H b/Source/WarpX.H index c2f5b94ac..8edf66f04 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -23,6 +23,7 @@ #include <PML.H> #include <BoostedFrameDiagnostic.H> #include <BilinearFilter.H> +#include <NCIGodfreyFilter.H> #ifdef WARPX_USE_PSATD #include <fftw3.h> @@ -99,6 +100,7 @@ public: // Back transformation diagnostic static bool do_boosted_frame_diagnostic; + static std::string lab_data_directory; static int num_snapshots_lab; static amrex::Real dt_snapshots_lab; static bool do_boosted_frame_fields; @@ -108,6 +110,8 @@ public: static amrex::Real gamma_boost; static amrex::Real beta_boost; static amrex::Vector<int> boost_direction; + static amrex::Real zmax_plasma_to_compute_max_step; + static int do_compute_max_step_from_zmax; static bool do_dynamic_scheduling; static bool refine_plasma; @@ -143,7 +147,9 @@ public: static amrex::IntVect filter_npass_each_dir; BilinearFilter bilinear_filter; - + amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_exeybz; + amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_bxbyez; + static int num_mirrors; amrex::Vector<amrex::Real> mirror_z; amrex::Vector<amrex::Real> mirror_z_width; @@ -152,9 +158,13 @@ public: void applyMirrors(amrex::Real time); void ComputeDt (); + // Compute max_step automatically for simulations in a boosted frame. + void computeMaxStepBoostAccelerator(amrex::Geometry geom); int MoveWindow (bool move_j); void UpdatePlasmaInjectionPosition (amrex::Real dt); + void ResetProbDomain (const amrex::RealBox& rb); + void EvolveE ( amrex::Real dt); void EvolveE (int lev, amrex::Real dt); void EvolveB ( amrex::Real dt); @@ -235,6 +245,12 @@ public: static int do_moving_window; static int moving_window_dir; + // slice generation // + void InitializeSliceMultiFabs (); + void SliceGenerationForDiagnostics (); + void WriteSlicePlotFile () const; + void ClearSliceMultiFabs (); + protected: //! Tagging cells for refinement @@ -476,10 +492,10 @@ private: amrex::Real current_injection_position = 0; // Plasma injection parameters - int do_plasma_injection = 0; + int warpx_do_continuous_injection = 0; int num_injected_species = -1; amrex::Vector<int> injected_plasma_species; - + int do_electrostatic = 0; int n_buffer = 4; amrex::Real const_dt = 0.5e-11; @@ -515,6 +531,9 @@ private: bool dump_plotfiles = true; bool dump_openpmd = false; #endif + bool plot_E_field = true; + bool plot_B_field = true; + bool plot_J_field = true; bool plot_part_per_cell = true; bool plot_part_per_grid = false; bool plot_part_per_proc = false; @@ -543,8 +562,16 @@ private: bool is_synchronized = true; - amrex::Vector<std::string> particle_plot_vars; - amrex::Vector<int> particle_plot_flags; + //Slice Parameters + int slice_max_grid_size; + int slice_plot_int = -1; + amrex::RealBox slice_realbox; + amrex::IntVect slice_cr_ratio; + amrex::Vector< std::unique_ptr<amrex::MultiFab> > F_slice; + amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice; #ifdef WARPX_USE_PSATD // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 27642c850..65079428b 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -17,6 +17,7 @@ #include <WarpXConst.H> #include <WarpXWrappers.h> #include <WarpXUtil.H> +#include <WarpXAlgorithmSelection.H> #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> @@ -32,12 +33,14 @@ int WarpX::moving_window_dir = -1; Real WarpX::gamma_boost = 1.; Real WarpX::beta_boost = 0.; Vector<int> WarpX::boost_direction = {0,0,0}; +int WarpX::do_compute_max_step_from_zmax = 0; +Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::current_deposition_algo = 3; -long WarpX::charge_deposition_algo = 0; -long WarpX::field_gathering_algo = 1; -long WarpX::particle_pusher_algo = 0; -int WarpX::maxwell_fdtd_solver_id = 0; +long WarpX::current_deposition_algo; +long WarpX::charge_deposition_algo; +long WarpX::field_gathering_algo; +long WarpX::particle_pusher_algo; +int WarpX::maxwell_fdtd_solver_id; long WarpX::nox = 1; long WarpX::noy = 1; @@ -55,6 +58,7 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; bool WarpX::do_boosted_frame_diagnostic = false; +std::string WarpX::lab_data_directory = "lab_frame_data"; int WarpX::num_snapshots_lab = std::numeric_limits<int>::lowest(); Real WarpX::dt_snapshots_lab = std::numeric_limits<Real>::lowest(); bool WarpX::do_boosted_frame_fields = true; @@ -145,15 +149,17 @@ WarpX::WarpX () // Particle Container mypc = std::unique_ptr<MultiParticleContainer> (new MultiParticleContainer(this)); - - if (do_plasma_injection) { - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast<PhysicalParticleContainer&>(pc); - ppc.injected = true; + warpx_do_continuous_injection = mypc->doContinuousInjection(); + if (warpx_do_continuous_injection){ + if (moving_window_v >= 0){ + // Inject particles continuously from the right end of the box + current_injection_position = geom[0].ProbHi(moving_window_dir); + } else { + // Inject particles continuously from the left end of the box + current_injection_position = geom[0].ProbLo(moving_window_dir); } } + do_boosted_frame_particles = mypc->doBoostedFrameDiags(); Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); @@ -223,6 +229,11 @@ WarpX::WarpX () #ifdef BL_USE_SENSEI_INSITU insitu_bridge = nullptr; #endif + + // NCI Godfrey filters can have different stencils + // at different levels (the stencil depends on c*dt/dz) + nci_godfrey_filter_exeybz.resize(nlevs_max); + nci_godfrey_filter_bxbyez.resize(nlevs_max); } WarpX::~WarpX () @@ -271,6 +282,12 @@ WarpX::ReadParameters () ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); + // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is + // specified by the user, 0 otherwise. + do_compute_max_step_from_zmax = + pp.query("zmax_plasma_to_compute_max_step", + zmax_plasma_to_compute_max_step); + pp.queryarr("B_external", B_external); pp.query("do_moving_window", do_moving_window); @@ -300,27 +317,14 @@ WarpX::ReadParameters () moving_window_v *= PhysConst::c; } - pp.query("do_plasma_injection", do_plasma_injection); - if (do_plasma_injection) { - pp.get("num_injected_species", num_injected_species); - injected_plasma_species.resize(num_injected_species); - pp.getarr("injected_plasma_species", injected_plasma_species, - 0, num_injected_species); - if (moving_window_v >= 0){ - // Inject particles continuously from the right end of the box - current_injection_position = geom[0].ProbHi(moving_window_dir); - } else { - // Inject particles continuously from the left end of the box - current_injection_position = geom[0].ProbLo(moving_window_dir); - } - } - pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); if (do_boosted_frame_diagnostic) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boosted frame diagnostic."); + pp.query("lab_data_directory", lab_data_directory); + std::string s; pp.get("boost_direction", s); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), @@ -331,8 +335,6 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); pp.query("do_boosted_frame_fields", do_boosted_frame_fields); - pp.query("do_boosted_frame_particles", do_boosted_frame_particles); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); @@ -385,6 +387,10 @@ WarpX::ReadParameters () if (ParallelDescriptor::NProcs() == 1) { plot_proc_number = false; } + // Fields to dump into plotfiles + pp.query("plot_E_field" , plot_E_field); + pp.query("plot_B_field" , plot_B_field); + pp.query("plot_J_field" , plot_J_field); pp.query("plot_part_per_cell", plot_part_per_cell); pp.query("plot_part_per_grid", plot_part_per_grid); pp.query("plot_part_per_proc", plot_part_per_proc); @@ -394,6 +400,7 @@ WarpX::ReadParameters () pp.query("plot_rho" , plot_rho); pp.query("plot_F" , plot_F); pp.query("plot_coarsening_ratio", plot_coarsening_ratio); + // Check that the coarsening_ratio can divide the blocking factor for (int lev=0; lev<maxLevel(); lev++){ for (int comp=0; comp<AMREX_SPACEDIM; comp++){ @@ -508,29 +515,12 @@ WarpX::ReadParameters () } { - ParmParse pp("algo"); - pp.query("current_deposition", current_deposition_algo); - pp.query("charge_deposition", charge_deposition_algo); - pp.query("field_gathering", field_gathering_algo); - pp.query("particle_pusher", particle_pusher_algo); - std::string s_solver = ""; - pp.query("maxwell_fdtd_solver", s_solver); - std::transform(s_solver.begin(), - s_solver.end(), - s_solver.begin(), - ::tolower); - // if maxwell_fdtd_solver is specified, set the value - // of maxwell_fdtd_solver_id accordingly. - // Otherwise keep the default value maxwell_fdtd_solver_id=0 - if (s_solver != "") { - if (s_solver == "yee") { - maxwell_fdtd_solver_id = 0; - } else if (s_solver == "ckc") { - maxwell_fdtd_solver_id = 1; - } else { - amrex::Abort("Unknown FDTD Solver type " + s_solver); - } - } + ParmParse pp("algo"); + current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); + charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); + field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); + particle_pusher_algo = GetAlgorithmInteger(pp, "particle_pusher"); + maxwell_fdtd_solver_id = GetAlgorithmInteger(pp, "maxwell_fdtd_solver"); } #ifdef WARPX_USE_PSATD @@ -559,6 +549,34 @@ WarpX::ReadParameters () pp.query("config", insitu_config); pp.query("pin_mesh", insitu_pin_mesh); } + + // for slice generation // + { + ParmParse pp("slice"); + amrex::Vector<Real> slice_lo(AMREX_SPACEDIM); + amrex::Vector<Real> slice_hi(AMREX_SPACEDIM); + Vector<int> slice_crse_ratio(AMREX_SPACEDIM); + // set default slice_crse_ratio // + for (int idim=0; idim < AMREX_SPACEDIM; ++idim ) + { + slice_crse_ratio[idim] = 1; + } + pp.queryarr("dom_lo",slice_lo,0,AMREX_SPACEDIM); + pp.queryarr("dom_hi",slice_hi,0,AMREX_SPACEDIM); + pp.queryarr("coarsening_ratio",slice_crse_ratio,0,AMREX_SPACEDIM); + pp.query("plot_int",slice_plot_int); + slice_realbox.setLo(slice_lo); + slice_realbox.setHi(slice_hi); + slice_cr_ratio = IntVect(AMREX_D_DECL(1,1,1)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + if (slice_crse_ratio[idim] > 1 ) { + slice_cr_ratio[idim] = slice_crse_ratio[idim]; + } + } + + } + } // This is a virtual function. @@ -662,7 +680,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number int ngz; if (WarpX::use_fdtd_nci_corr) { - int ng = ngz_tmp + (mypc->nstencilz_fdtd_nci_corr-1); + int ng = ngz_tmp + NCIGodfreyFilter::stencil_width; ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; @@ -913,8 +931,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (load_balance_int > 0) { costs[lev].reset(new MultiFab(ba, dm, 1, 0)); } - - } std::array<Real,3> @@ -1035,12 +1051,19 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, for (MFIter mfi(divE, true); mfi.isValid(); ++mfi) { const Box& bx = mfi.tilebox(); +#ifdef WARPX_RZ + const Real xmin = GetInstance().Geom(0).ProbLo(0); +#endif WRPX_COMPUTE_DIVE(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_N_ANYD(divE[mfi],dcomp), BL_TO_FORTRAN_ANYD((*E[0])[mfi]), BL_TO_FORTRAN_ANYD((*E[1])[mfi]), BL_TO_FORTRAN_ANYD((*E[2])[mfi]), - dx.data()); + dx.data() +#ifdef WARPX_RZ + ,&xmin +#endif + ); } } @@ -1055,19 +1078,25 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, for (MFIter mfi(divE, true); mfi.isValid(); ++mfi) { Box bx = mfi.growntilebox(ngrow); +#ifdef WARPX_RZ + const Real xmin = GetInstance().Geom(0).ProbLo(0); +#endif WRPX_COMPUTE_DIVE(bx.loVect(), bx.hiVect(), BL_TO_FORTRAN_N_ANYD(divE[mfi],dcomp), BL_TO_FORTRAN_ANYD((*E[0])[mfi]), BL_TO_FORTRAN_ANYD((*E[1])[mfi]), BL_TO_FORTRAN_ANYD((*E[2])[mfi]), - dx.data()); + dx.data() +#ifdef WARPX_RZ + ,&xmin +#endif + ); } } void WarpX::BuildBufferMasks () { - int ngbuffer = std::max(n_field_gather_buffer, n_current_deposition_buffer); for (int lev = 1; lev <= maxLevel(); ++lev) { for (int ipass = 0; ipass < 2; ++ipass) |