diff options
Diffstat (limited to 'Source')
24 files changed, 980 insertions, 300 deletions
diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BackTransformedDiagnostic.H index 5d95aaf7d..9e24caa1b 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BackTransformedDiagnostic.H @@ -1,5 +1,5 @@ -#ifndef WARPX_BoostedFrameDiagnostic_H_ -#define WARPX_BoostedFrameDiagnostic_H_ +#ifndef WARPX_BackTransformedDiagnostic_H_ +#define WARPX_BackTransformedDiagnostic_H_ #include <vector> #include <map> @@ -150,7 +150,7 @@ class LabFrameSlice : public LabFrameDiag { }; /** \brief - * BoostedFrameDiagnostic class handles the back-transformation of data when + * BackTransformedDiagnostic class handles the back-transformation of data when * running simulations in a boosted frame of reference to the lab-frame. * Because of the relativity of simultaneity, events that are synchronized * in the simulation boosted frame are not @@ -163,13 +163,13 @@ class LabFrameSlice : public LabFrameDiag { * to the output directory. The functions Flush() and writeLabFrameData() * are called at the end of the simulation and when the * the buffer for data storage is full, respectively. The particle data - * is collected and written only if particle.do_boosted_frame_diagnostic = 1. + * is collected and written only if particle.do_back_transformed_diagnostics = 1. */ -class BoostedFrameDiagnostic { +class BackTransformedDiagnostic { public: - BoostedFrameDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab, + BackTransformedDiagnostic(amrex::Real zmin_lab, amrex::Real zmax_lab, amrex::Real v_window_lab, amrex::Real dt_snapshots_lab, int N_snapshots, amrex::Real dt_slice_snapshots_lab, int N_slice_snapshots, amrex::Real gamma_boost, diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BackTransformedDiagnostic.cpp index 297b4f5be..2880b37b1 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BackTransformedDiagnostic.cpp @@ -1,7 +1,7 @@ #include <AMReX_MultiFabUtil.H> #include <AMReX_MultiFabUtil_C.H> -#include "BoostedFrameDiagnostic.H" +#include "BackTransformedDiagnostic.H" #include "SliceDiagnostic.H" #include "WarpX_f.H" #include "WarpX.H" @@ -514,8 +514,8 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) } } -BoostedFrameDiagnostic:: -BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, +BackTransformedDiagnostic:: +BackTransformedDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, Real dt_snapshots_lab, int N_snapshots, Real dt_slice_snapshots_lab, int N_slice_snapshots, Real gamma_boost, Real t_boost, Real dt_boost, @@ -531,10 +531,10 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, { - BL_PROFILE("BoostedFrameDiagnostic::BoostedFrameDiagnostic"); + BL_PROFILE("BackTransformedDiagnostic::BackTransformedDiagnostic"); - AMREX_ALWAYS_ASSERT(WarpX::do_boosted_frame_fields or - WarpX::do_boosted_frame_particles); + AMREX_ALWAYS_ASSERT(WarpX::do_back_transformed_fields or + WarpX::do_back_transformed_particles); inv_gamma_boost_ = 1.0 / gamma_boost_; beta_boost_ = std::sqrt(1.0 - inv_gamma_boost_*inv_gamma_boost_); @@ -679,9 +679,9 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, AMREX_ALWAYS_ASSERT(max_box_size_ >= num_buffer_); } -void BoostedFrameDiagnostic::Flush(const Geometry& geom) +void BackTransformedDiagnostic::Flush(const Geometry& geom) { - BL_PROFILE("BoostedFrameDiagnostic::Flush"); + BL_PROFILE("BackTransformedDiagnostic::Flush"); VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -696,7 +696,7 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) int i_lab = (LabFrameDiags_[i]->current_z_lab - zmin_lab) / dz_lab_; if (LabFrameDiags_[i]->buff_counter_ != 0) { - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { const BoxArray& ba = LabFrameDiags_[i]->data_buffer_->boxArray(); const int hi = ba[0].bigEnd(boost_direction_); const int lo = hi - LabFrameDiags_[i]->buff_counter_ + 1; @@ -731,12 +731,12 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) #endif } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_back_transformed_particles) { // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) { // Get species name std::string species_name = - species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; + species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)]; #ifdef WARPX_USE_HDF5 // Dump species data writeParticleDataHDF5(LabFrameDiags_[i]->particles_buffer_[j], @@ -765,12 +765,12 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeLabFrameData(const MultiFab* cell_centered_data, const MultiParticleContainer& mypc, const Geometry& geom, const Real t_boost, const Real dt) { - BL_PROFILE("BoostedFrameDiagnostic::writeLabFrameData"); + BL_PROFILE("BackTransformedDiagnostic::writeLabFrameData"); VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -808,7 +808,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // If buffer of snapshot i is empty... if ( LabFrameDiags_[i]->buff_counter_ == 0) { // ... reset fields buffer data_buffer_ - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { LabFrameDiags_[i]->buff_box_.setSmall(boost_direction_, i_lab - num_buffer_ + 1); LabFrameDiags_[i]->buff_box_.setBig(boost_direction_, i_lab); @@ -820,12 +820,12 @@ writeLabFrameData(const MultiFab* cell_centered_data, buff_dm, ncomp_to_dump, 0) ); } // ... reset particle buffer particles_buffer_[i] - if (WarpX::do_boosted_frame_particles) + if (WarpX::do_back_transformed_particles) LabFrameDiags_[i]->particles_buffer_.resize( - mypc.nSpeciesBoostedFrameDiags()); + mypc.nSpeciesBackTransformedDiagnostics()); } - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { const int ncomp = cell_centered_data->nComp(); const int start_comp = 0; const bool interpolate = true; @@ -873,7 +873,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, tmp_slice_ptr.reset(nullptr); } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_back_transformed_particles) { if (LabFrameDiags_[i]->t_lab != prev_t_lab ) { if (tmp_particle_buffer.size()>0) @@ -881,7 +881,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, tmp_particle_buffer.clear(); tmp_particle_buffer.shrink_to_fit(); } - tmp_particle_buffer.resize(mypc.nSpeciesBoostedFrameDiags()); + tmp_particle_buffer.resize(mypc.nSpeciesBackTransformedDiagnostics()); mypc.GetLabFrameData( LabFrameDiags_[i]->file_name, i_lab, boost_direction_, old_z_boost, LabFrameDiags_[i]->current_z_boost, @@ -889,7 +889,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, tmp_particle_buffer); } LabFrameDiags_[i]->AddPartDataToParticleBuffer(tmp_particle_buffer, - mypc.nSpeciesBoostedFrameDiags()); + mypc.nSpeciesBackTransformedDiagnostics()); } ++LabFrameDiags_[i]->buff_counter_; @@ -898,7 +898,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, // If buffer full, write to disk. if ( LabFrameDiags_[i]->buff_counter_ == num_buffer_) { - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { #ifdef WARPX_USE_HDF5 Box buff_box = LabFrameDiags_[i]->buff_box_; @@ -916,12 +916,12 @@ writeLabFrameData(const MultiFab* cell_centered_data, #endif } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_back_transformed_particles) { // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) { + for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) { // Get species name const std::string species_name = species_names[ - mypc.mapSpeciesBoostedFrameDiags(j)]; + mypc.mapSpeciesBackTransformedDiagnostics(j)]; #ifdef WARPX_USE_HDF5 // Write data to disk (HDF5) writeParticleDataHDF5(LabFrameDiags_[i]->particles_buffer_[j], @@ -949,7 +949,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, #ifdef WARPX_USE_HDF5 void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const std::string& species_name) { @@ -997,11 +997,11 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat #endif void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, const std::string& name, const int i_lab) { - BL_PROFILE("BoostedFrameDiagnostic::writeParticleData"); + BL_PROFILE("BackTransformedDiagnostic::writeParticleData"); std::string field_name; std::ofstream ofs; @@ -1047,10 +1047,10 @@ writeParticleData(const WarpXParticleContainer::DiagnosticParticleData& pdata, } void -BoostedFrameDiagnostic:: +BackTransformedDiagnostic:: writeMetaData () { - BL_PROFILE("BoostedFrameDiagnostic::writeMetaData"); + BL_PROFILE("BackTransformedDiagnostic::writeMetaData"); if (ParallelDescriptor::IOProcessor()) { const std::string fullpath = WarpX::lab_data_directory + "/snapshots"; @@ -1134,7 +1134,7 @@ LabFrameSnapShot(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in, file_num, 5); createLabFrameDirectories(); buff_counter_ = 0; - if (WarpX::do_boosted_frame_fields) data_buffer_.reset(nullptr); + if (WarpX::do_back_transformed_fields) data_buffer_.reset(nullptr); } void @@ -1158,7 +1158,7 @@ createLabFrameDirectories() { if (ParallelDescriptor::IOProcessor()) { - if (WarpX::do_boosted_frame_fields) + if (WarpX::do_back_transformed_fields) { const auto lo = lbound(buff_box_); for (int comp = 0; comp < ncomp_to_dump_; ++comp) { @@ -1176,15 +1176,15 @@ createLabFrameDirectories() { ParallelDescriptor::Barrier(); - if (WarpX::do_boosted_frame_particles){ + if (WarpX::do_back_transformed_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector<std::string> species_names = mypc.GetSpeciesNames(); // Loop over species to be dumped to BFD - for (int j = 0; j < mypc.nSpeciesBoostedFrameDiags(); ++j) + for (int j = 0; j < mypc.nSpeciesBackTransformedDiagnostics(); ++j) { // Loop over species to be dumped to BFD std::string species_name = - species_names[mypc.mapSpeciesBoostedFrameDiags(j)]; + species_names[mypc.mapSpeciesBackTransformedDiagnostics(j)]; output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast<int>(particle_field_names.size()); ++k) { @@ -1211,10 +1211,10 @@ createLabFrameDirectories() { const std::string particles_prefix = "particle"; // Loop over species to be dumped to BFD - for(int i = 0; i < mypc.nSpeciesBoostedFrameDiags(); ++i) { + for(int i = 0; i < mypc.nSpeciesBackTransformedDiagnostics(); ++i) { // Get species name std::string species_name = - species_names[mypc.mapSpeciesBoostedFrameDiags(i)]; + species_names[mypc.mapSpeciesBackTransformedDiagnostics(i)]; const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); @@ -1302,7 +1302,7 @@ LabFrameSlice(Real t_lab_in, Real t_boost, Real inv_gamma_boost_in, dx_ = cell_dx; dy_ = cell_dy; - if (WarpX::do_boosted_frame_fields) data_buffer_.reset(nullptr); + if (WarpX::do_back_transformed_fields) data_buffer_.reset(nullptr); } void diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package index dfd947d53..b624d6ebe 100644 --- a/Source/Diagnostics/Make.package +++ b/Source/Diagnostics/Make.package @@ -1,9 +1,9 @@ CEXE_sources += WarpXIO.cpp -CEXE_sources += BoostedFrameDiagnostic.cpp +CEXE_sources += BackTransformedDiagnostic.cpp CEXE_sources += ParticleIO.cpp CEXE_sources += FieldIO.cpp CEXE_headers += FieldIO.H -CEXE_headers += BoostedFrameDiagnostic.H +CEXE_headers += BackTransformedDiagnostic.H CEXE_headers += ElectrostaticIO.cpp CEXE_headers += SliceDiagnostic.H CEXE_sources += SliceDiagnostic.cpp diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 7a3262703..b5fd52bdc 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -141,9 +141,9 @@ WarpX::EvolveEM (int numsteps) bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { std::unique_ptr<MultiFab> cell_centered_data = nullptr; - if (WarpX::do_boosted_frame_fields) { + if (WarpX::do_back_transformed_fields) { cell_centered_data = GetCellCenteredData(); } myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); @@ -261,7 +261,7 @@ WarpX::EvolveEM (int numsteps) WriteCheckPointFile(); } - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { myBFD->Flush(geom[0]); } diff --git a/Source/Initialization/InjectorPosition.H b/Source/Initialization/InjectorPosition.H index 6ecae93e0..4ab2fa022 100644 --- a/Source/Initialization/InjectorPosition.H +++ b/Source/Initialization/InjectorPosition.H @@ -29,21 +29,25 @@ struct InjectorPositionRegular // is a_ppc*(ref_fac**AMREX_SPACEDIM). AMREX_GPU_HOST_DEVICE amrex::XDim3 - getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + getPositionUnitBox (int const i_part, int const ref_fac=1) const noexcept { - int nx = ref_fac*ppc.x; - int ny = ref_fac*ppc.y; + using namespace amrex; + + int const nx = ref_fac*ppc.x; + int const ny = ref_fac*ppc.y; #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ) - int nz = ref_fac*ppc.z; + int const nz = ref_fac*ppc.z; #else - int nz = 1; + int const nz = 1; #endif - int ix_part = i_part/(ny*nz); // written this way backward compatibility - int iz_part = (i_part-ix_part*(ny*nz)) / ny; - int iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; - return amrex::XDim3{(amrex::Real(0.5)+ix_part)/nx, - (amrex::Real(0.5)+iy_part)/ny, - (amrex::Real(0.5)+iz_part) / nz}; + int const ix_part = i_part / (ny*nz); // written this way backward compatibility + int const iz_part = (i_part-ix_part*(ny*nz)) / ny; + int const iy_part = (i_part-ix_part*(ny*nz)) - ny*iz_part; + return XDim3{ + (0.5_rt + ix_part) / nx, + (0.5_rt + iy_part) / ny, + (0.5_rt + iz_part) / nz + }; } private: amrex::Dim3 ppc; @@ -100,7 +104,7 @@ struct InjectorPosition // (the union is called Object, and the instance is called object). AMREX_GPU_HOST_DEVICE amrex::XDim3 - getPositionUnitBox (int i_part, int ref_fac=1) const noexcept + getPositionUnitBox (int const i_part, int const ref_fac=1) const noexcept { switch (type) { diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index a2e59c177..0814f369b 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -82,14 +82,14 @@ WarpX::InitData () void WarpX::InitDiagnostics () { - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { const Real* current_lo = geom[0].ProbLo(); const Real* current_hi = geom[0].ProbHi(); Real dt_boost = dt[0]; // Find the positions of the lab-frame box that corresponds to the boosted-frame box at t=0 Real zmin_lab = current_lo[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); Real zmax_lab = current_hi[moving_window_dir]/( (1.+beta_boost)*gamma_boost ); - myBFD.reset(new BoostedFrameDiagnostic(zmin_lab, + myBFD.reset(new BackTransformedDiagnostic(zmin_lab, zmax_lab, moving_window_v, dt_snapshots_lab, num_snapshots_lab, diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 8571c74ad..9493672e0 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -26,7 +26,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, { charge = 1.0; mass = std::numeric_limits<Real>::max(); - do_boosted_frame_diags = 0; + do_back_transformed_diagnostics = 0; ParmParse pp(laser_name); @@ -100,7 +100,7 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, } // Plane normal - Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); + Real s = 1.0_rt / std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]); nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s }; if (WarpX::gamma_boost > 1.) { @@ -119,19 +119,19 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, } // The first polarization vector - s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]); + s = 1.0_rt / std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]); p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s }; - Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); + Real const dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, "Laser plane vector is not perpendicular to the main polarization vector"); p_Y = CrossProduct(nvec, p_X); // The second polarization vector - s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); + s = 1.0_rt / std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]); stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s }; - dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14, + Real const dp2 = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp2) < 1.0e-14, "stc_direction is not perpendicular to the laser plane vector"); // Get angle between p_X and stc_direction @@ -266,20 +266,20 @@ LaserParticleContainer::InitData (int lev) position = updated_position; } - auto Transform = [&](int i, int j) -> Vector<Real>{ + auto Transform = [&](int const i, int const 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], - position[1] + (S_X*(i+0.5))*u_X[1] + (S_Y*(j+0.5))*u_Y[1], - position[2] + (S_X*(i+0.5))*u_X[2] + (S_Y*(j+0.5))*u_Y[2] }; + return { position[0] + (S_X*(Real(i)+0.5_rt))*u_X[0] + (S_Y*(Real(j)+0.5_rt))*u_Y[0], + position[1] + (S_X*(Real(i)+0.5_rt))*u_X[1] + (S_Y*(Real(j)+0.5_rt))*u_Y[1], + position[2] + (S_X*(Real(i)+0.5_rt))*u_X[2] + (S_Y*(Real(j)+0.5_rt))*u_Y[2] }; #else # if (defined WARPX_DIM_RZ) - return { position[0] + (S_X*(i+0.5)), + return { position[0] + (S_X*(Real(i)+0.5)), 0.0, position[2]}; # else - return { position[0] + (S_X*(i+0.5))*u_X[0], + return { position[0] + (S_X*(Real(i)+0.5))*u_X[0], 0.0, - position[2] + (S_X*(i+0.5))*u_X[2] }; + position[2] + (S_X*(Real(i)+0.5))*u_X[2] }; # endif #endif }; @@ -449,9 +449,9 @@ LaserParticleContainer::Evolve (int lev, #endif { #ifdef _OPENMP - int thread_num = omp_get_thread_num(); + int const thread_num = omp_get_thread_num(); #else - int thread_num = 0; + int const thread_num = 0; #endif Cuda::ManagedDeviceVector<Real> plane_Xp, plane_Yp, amplitude_E; @@ -610,7 +610,7 @@ void LaserParticleContainer::ComputeWeightMobility (Real Sx, Real Sy) { constexpr Real eps = 0.01; - constexpr Real fac = 1.0/(2.0*MathConst::pi*PhysConst::mu0*PhysConst::c*PhysConst::c*eps); + constexpr Real fac = 1.0_rt / (2.0_rt * MathConst::pi * PhysConst::mu0 * PhysConst::c * PhysConst::c * eps); weight = fac * wavelength * Sx * Sy / std::min(Sx,Sy) * e_max; // The mobility is the constant of proportionality between the field to diff --git a/Source/Laser/LaserProfiles.cpp b/Source/Laser/LaserProfiles.cpp index 281ab2101..44411cedf 100644 --- a/Source/Laser/LaserProfiles.cpp +++ b/Source/Laser/LaserProfiles.cpp @@ -28,16 +28,16 @@ LaserParticleContainer::gaussian_laser_profile ( const Real oscillation_phase = k0 * PhysConst::c * ( t - profile_t_peak ); // The coefficients below contain info about Gouy phase, // laser diffraction, and phase front curvature - const Complex diffract_factor = Real(1.) + I * profile_focal_distance - * Real(2.)/( k0 * profile_waist * profile_waist ); - const Complex inv_complex_waist_2 = Real(1.)/( profile_waist*profile_waist * diffract_factor ); + const Complex diffract_factor = 1._rt + I * profile_focal_distance + * 2._rt/( k0 * profile_waist * profile_waist ); + const Complex inv_complex_waist_2 = 1._rt / ( profile_waist*profile_waist * diffract_factor ); // Time stretching due to STCs and phi2 complex envelope // (1 if zeta=0, beta=0, phi2=0) - const Complex stretch_factor = Real(1.) + Real(4.) * + const Complex stretch_factor = 1._rt + 4._rt * (zeta+beta*profile_focal_distance) * (zeta+beta*profile_focal_distance) * (inv_tau2*inv_complex_waist_2) + - Real(2.)*I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; + 2._rt *I*(phi2 - beta*beta*k0*profile_focal_distance) * inv_tau2; // Amplitude and monochromatic oscillations Complex prefactor = e_max * MathFunc::exp( I * oscillation_phase ); @@ -61,10 +61,10 @@ LaserParticleContainer::gaussian_laser_profile ( amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (int i) { - const Complex stc_exponent = Real(1.)/stretch_factor * inv_tau2 * + const Complex stc_exponent = 1._rt / stretch_factor * inv_tau2 * MathFunc::pow((t - tmp_profile_t_peak - tmp_beta*k0*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) - - Real(2.)*I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) + 2._rt *I*(Xp[i]*std::cos(tmp_theta_stc) + Yp[i]*std::sin(tmp_theta_stc)) *( tmp_zeta - tmp_beta*tmp_profile_focal_distance ) * inv_complex_waist_2),2); // stcfactor = everything but complex transverse envelope const Complex stcfactor = prefactor * MathFunc::exp( - stc_exponent ); diff --git a/Source/Parallelization/InterpolateCurrentFineToCoarse.H b/Source/Parallelization/InterpolateCurrentFineToCoarse.H index 148b725d0..cbbcdfab5 100644 --- a/Source/Parallelization/InterpolateCurrentFineToCoarse.H +++ b/Source/Parallelization/InterpolateCurrentFineToCoarse.H @@ -56,6 +56,8 @@ public: int const k ) const noexcept { + using namespace amrex; + auto const & fine_unsafe = m_fine; // out-of-bounds access not secured with zero-values yet auto const & coarse = m_coarse; // out-of-bounds access not secured but will also not occur @@ -71,29 +73,29 @@ public: int const kk = k * m_refinement_ratio; #if AMREX_SPACEDIM == 2 if (IDim == 0) { - coarse(i, j, k) = 0.25 * ( + coarse(i, j, k) = 0.25_rt * ( fine(ii, jj, kk) + fine(ii + 1, jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + fine(ii, jj + 1, kk) + fine(ii + 1, jj + 1, kk) ) ); } else if (IDim == 2) { - coarse(i, j, k) = 0.25 * ( + coarse(i, j, k) = 0.25_rt * ( fine(ii, jj, kk) + fine(ii, jj + 1, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii - 1, jj, kk) + fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj, kk) + fine(ii + 1, jj + 1, kk) ) ); } else { - coarse(i, j, k) = 0.25 * ( + coarse(i, j, k) = 0.25_rt * ( fine(ii, jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii - 1, jj , kk) + fine(ii + 1, jj , kk) + fine(ii , jj - 1, kk) + fine(ii , jj + 1, kk) ) + - 0.25 * ( + 0.25_rt * ( fine(ii - 1, jj - 1, kk) + fine(ii + 1, jj - 1, kk) + fine(ii - 1, jj + 1, kk) + fine(ii + 1, jj + 1, kk) ) @@ -101,64 +103,64 @@ public: } #elif AMREX_SPACEDIM == 3 if (IDim == 0) { - coarse(i,j,k) = 0.125 * ( + coarse(i,j,k) = 0.125_rt * ( fine(ii , jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) + fine(ii , jj , kk-1) + fine(ii , jj , kk+1) ) + - 0.25 * ( + 0.25_rt * ( fine(ii , jj-1, kk-1) + fine(ii , jj+1, kk-1) + fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) ) + fine(ii+1, jj, kk) + - 0.5 * ( + 0.5_rt * ( fine(ii+1, jj-1, kk ) + fine(ii+1, jj+1, kk ) + fine(ii+1, jj , kk-1) + fine(ii+1, jj , kk+1) ) + - 0.25 * ( + 0.25_rt * ( fine(ii+1, jj-1, kk-1) + fine(ii+1, jj+1, kk-1) + fine(ii+1, jj-1, kk+1) + fine(ii+1, jj+1, kk+1) ) ); } else if (IDim == 1) { - coarse(i, j, k) = 0.125 * ( + coarse(i, j, k) = 0.125_rt * ( fine(ii, jj , kk) + - 0.5 * ( + 0.5_rt * ( fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + fine(ii , jj , kk-1) + fine(ii , jj , kk+1) - ) + - 0.25 * ( + ) + + 0.25_rt * ( fine(ii-1, jj , kk-1) + fine(ii+1, jj , kk-1) + fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) - ) + - fine(ii, jj+1, kk) + - 0.5 * ( + ) + + fine(ii, jj+1, kk) + + 0.5_rt * ( fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) + fine(ii , jj+1, kk-1) + fine(ii , jj+1, kk+1) - ) + - 0.25 * ( + ) + + 0.25_rt * ( fine(ii-1, jj+1, kk-1) + fine(ii+1, jj+1, kk-1) + fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) - ) + ) ); } else { - coarse(i, j, k) = 0.125 * ( + coarse(i, j, k) = 0.125_rt * ( fine(ii, jj, kk ) + - 0.5 * ( + 0.5_rt * ( fine(ii-1, jj , kk ) + fine(ii+1, jj , kk ) + fine(ii , jj-1, kk ) + fine(ii , jj+1, kk ) ) + - 0.25 * ( + 0.25_rt * ( fine(ii-1, jj-1, kk ) + fine(ii+1, jj-1, kk ) + fine(ii-1, jj+1, kk ) + fine(ii+1, jj+1, kk ) ) + fine(ii, jj, kk+1) + - 0.5 * ( + 0.5_rt * ( fine(ii-1, jj , kk+1) + fine(ii+1, jj , kk+1) + fine(ii , jj-1, kk+1) + fine(ii , jj+1, kk+1) ) + - 0.25 * ( + 0.25_rt * ( fine(ii-1, jj-1, kk+1) + fine(ii+1, jj-1, kk+1) + fine(ii-1, jj+1, kk+1) + fine(ii+1, jj+1, kk+1) ) diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 92f0b4f09..52df3dc25 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -54,6 +54,157 @@ WarpX::UpdateAuxilaryData () { BL_PROFILE("UpdateAuxilaryData()"); + if (Bfield_aux[0][0]->ixType() == Bfield_fp[0][0]->ixType()) { + UpdateAuxilaryDataSameType(); + } else { + UpdateAuxilaryDataStagToNodal(); + } +} + +void +WarpX::UpdateAuxilaryDataStagToNodal () +{ + // For level 0, we only need to do the average. +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Bfield_aux[0][0]); mfi.isValid(); ++mfi) + { + Array4<Real> const& bx_aux = Bfield_aux[0][0]->array(mfi); + Array4<Real> const& by_aux = Bfield_aux[0][1]->array(mfi); + Array4<Real> const& bz_aux = Bfield_aux[0][2]->array(mfi); + Array4<Real const> const& bx_fp = Bfield_fp[0][0]->const_array(mfi); + Array4<Real const> const& by_fp = Bfield_fp[0][1]->const_array(mfi); + Array4<Real const> const& bz_fp = Bfield_fp[0][2]->const_array(mfi); + + Array4<Real> const& ex_aux = Efield_aux[0][0]->array(mfi); + Array4<Real> const& ey_aux = Efield_aux[0][1]->array(mfi); + Array4<Real> const& ez_aux = Efield_aux[0][2]->array(mfi); + Array4<Real const> const& ex_fp = Efield_fp[0][0]->const_array(mfi); + Array4<Real const> const& ey_fp = Efield_fp[0][1]->const_array(mfi); + Array4<Real const> const& ez_fp = Efield_fp[0][2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp); + warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp); + warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp); + warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp); + warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp); + warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp); + }); + } + + for (int lev = 1; lev <= finest_level; ++lev) + { + BoxArray const& nba = Bfield_aux[lev][0]->boxArray(); + BoxArray const& cnba = amrex::coarsen(nba,2); + DistributionMapping const& dm = Bfield_aux[lev][0]->DistributionMap(); + auto const& cperiod = Geom(lev-1).periodicity(); + + // Bfield + { + Array<std::unique_ptr<MultiFab>,3> Btmp; + if (Bfield_cax[lev][0]) { + for (int i = 0; i < 3; ++i) { + Btmp[i].reset(new MultiFab(*Bfield_cax[lev][i], amrex::make_alias, 0, 1)); + } + } else { + IntVect ngtmp = Bfield_aux[lev-1][0]->nGrowVect(); + for (int i = 0; i < 3; ++i) { + Btmp[i].reset(new MultiFab(cnba, dm, 1, ngtmp)); + } + } + // ParallelCopy from coarse level + for (int i = 0; i < 3; ++i) { + IntVect ng = Btmp[i]->nGrowVect(); + Btmp[i]->ParallelCopy(*Bfield_aux[lev-1][i], 0, 0, 1, ng, ng, cperiod); + } + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Bfield_aux[lev][0]); mfi.isValid(); ++mfi) + { + Array4<Real> const& bx_aux = Bfield_aux[lev][0]->array(mfi); + Array4<Real> const& by_aux = Bfield_aux[lev][1]->array(mfi); + Array4<Real> const& bz_aux = Bfield_aux[lev][2]->array(mfi); + Array4<Real const> const& bx_fp = Bfield_fp[lev][0]->const_array(mfi); + Array4<Real const> const& by_fp = Bfield_fp[lev][1]->const_array(mfi); + Array4<Real const> const& bz_fp = Bfield_fp[lev][2]->const_array(mfi); + Array4<Real const> const& bx_cp = Bfield_cp[lev][0]->const_array(mfi); + Array4<Real const> const& by_cp = Bfield_cp[lev][1]->const_array(mfi); + Array4<Real const> const& bz_cp = Bfield_cp[lev][2]->const_array(mfi); + Array4<Real const> const& bx_c = Btmp[0]->const_array(mfi); + Array4<Real const> const& by_c = Btmp[1]->const_array(mfi); + Array4<Real const> const& bz_c = Btmp[2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp_nd_bfield_x(j,k,l, bx_aux, bx_fp, bx_cp, bx_c); + warpx_interp_nd_bfield_y(j,k,l, by_aux, by_fp, by_cp, by_c); + warpx_interp_nd_bfield_z(j,k,l, bz_aux, bz_fp, bz_cp, bz_c); + }); + } + } + + // Efield + { + Array<std::unique_ptr<MultiFab>,3> Etmp; + if (Efield_cax[lev][0]) { + for (int i = 0; i < 3; ++i) { + Etmp[i].reset(new MultiFab(*Efield_cax[lev][i], amrex::make_alias, 0, 1)); + } + } else { + IntVect ngtmp = Efield_aux[lev-1][0]->nGrowVect(); + for (int i = 0; i < 3; ++i) { + Etmp[i].reset(new MultiFab(cnba, dm, 1, ngtmp)); + } + } + // ParallelCopy from coarse level + for (int i = 0; i < 3; ++i) { + IntVect ng = Etmp[i]->nGrowVect(); + Etmp[i]->ParallelCopy(*Efield_aux[lev-1][i], 0, 0, 1, ng, ng, cperiod); + } + +#ifdef _OPENMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(*Efield_aux[lev][0]); mfi.isValid(); ++mfi) + { + Array4<Real> const& ex_aux = Efield_aux[lev][0]->array(mfi); + Array4<Real> const& ey_aux = Efield_aux[lev][1]->array(mfi); + Array4<Real> const& ez_aux = Efield_aux[lev][2]->array(mfi); + Array4<Real const> const& ex_fp = Efield_fp[lev][0]->const_array(mfi); + Array4<Real const> const& ey_fp = Efield_fp[lev][1]->const_array(mfi); + Array4<Real const> const& ez_fp = Efield_fp[lev][2]->const_array(mfi); + Array4<Real const> const& ex_cp = Efield_cp[lev][0]->const_array(mfi); + Array4<Real const> const& ey_cp = Efield_cp[lev][1]->const_array(mfi); + Array4<Real const> const& ez_cp = Efield_cp[lev][2]->const_array(mfi); + Array4<Real const> const& ex_c = Etmp[0]->const_array(mfi); + Array4<Real const> const& ey_c = Etmp[1]->const_array(mfi); + Array4<Real const> const& ez_c = Etmp[2]->const_array(mfi); + + const Box& bx = mfi.fabbox(); + amrex::ParallelFor(bx, + [=] AMREX_GPU_DEVICE (int j, int k, int l) noexcept + { + warpx_interp_nd_efield_x(j,k,l, ex_aux, ex_fp, ex_cp, ex_c); + warpx_interp_nd_efield_y(j,k,l, ey_aux, ey_fp, ey_cp, ey_c); + warpx_interp_nd_efield_z(j,k,l, ez_aux, ez_fp, ez_cp, ez_c); + }); + } + } + } +} + +void +WarpX::UpdateAuxilaryDataSameType () +{ for (int lev = 1; lev <= finest_level; ++lev) { const auto& crse_period = Geom(lev-1).periodicity(); diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index 093323ec3..169cd0ee1 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -5,38 +5,38 @@ AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_bfield_x (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bxa, + amrex::Array4<amrex::Real > const& Bxa, amrex::Array4<amrex::Real const> const& Bxf, amrex::Array4<amrex::Real const> const& Bxc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wx = (j == jg*2) ? 0.0 : 0.5; - Real owx = 1.0-wx; + Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt; + Real const owx = 1.0_rt - wx; Bxa(j,k,l) = owx * Bxc(jg,kg,lg) + wx * Bxc(jg+1,kg,lg) + Bxf(j,k,l); } AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_bfield_y (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bya, + amrex::Array4<amrex::Real > const& Bya, amrex::Array4<amrex::Real const> const& Byf, amrex::Array4<amrex::Real const> const& Byc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); // Note that for 2d, l=0, because the amrex convention is used here. #if (AMREX_SPACEDIM == 3) - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; Bya(j,k,l) = owy * Byc(jg,kg,lg) + wy * Byc(jg,kg+1,lg) + Byf(j,k,l); #else Bya(j,k,l) = Byc(jg,kg,lg) + Byf(j,k,l); @@ -45,47 +45,47 @@ void warpx_interp_bfield_y (int j, int k, int l, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_bfield_z (int j, int k, int l, - amrex::Array4<amrex::Real> const& Bza, + amrex::Array4<amrex::Real > const& Bza, amrex::Array4<amrex::Real const> const& Bzf, amrex::Array4<amrex::Real const> const& Bzc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); // Note that for 2d, l=0, because the amrex convention is used here. #if (AMREX_SPACEDIM == 3) - Real wz = (l == lg*2) ? 0.0 : 0.5; - Real owz = 1.0-wz; + Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt; + Real const owz = 1.0_rt - wz; Bza(j,k,l) = owz * Bzc(jg,kg,lg) + owz * Bzc(jg,kg,lg+1) + Bzf(j,k,l); #else - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; Bza(j,k,l) = owy * Bzc(jg,kg,lg) + owy * Bzc(jg,kg+1,lg) + Bzf(j,k,l); #endif } AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_efield_x (int j, int k, int l, - amrex::Array4<amrex::Real> const& Exa, + amrex::Array4<amrex::Real > const& Exa, amrex::Array4<amrex::Real const> const& Exf, amrex::Array4<amrex::Real const> const& Exc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; #if (AMREX_SPACEDIM == 3) - Real wz = (l == lg*2) ? 0.0 : 0.5; - Real owz = 1.0-wz; + Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt; + Real const owz = 1.0_rt - wz; Exa(j,k,l) = owy * owz * Exc(jg ,kg ,lg ) + wy * owz * Exc(jg ,kg+1,lg ) + owy * wz * Exc(jg ,kg ,lg+1) @@ -98,30 +98,30 @@ void warpx_interp_efield_x (int j, int k, int l, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_efield_y (int j, int k, int l, - amrex::Array4<amrex::Real> const& Eya, + amrex::Array4<amrex::Real > const& Eya, amrex::Array4<amrex::Real const> const& Eyf, amrex::Array4<amrex::Real const> const& Eyc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wx = (j == jg*2) ? 0.0 : 0.5; - Real owx = 1.0-wx; + Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt; + Real const owx = 1.0_rt - wx; #if (AMREX_SPACEDIM == 3) - Real wz = (l == lg*2) ? 0.0 : 0.5; - Real owz = 1.0-wz; + Real const wz = (l == lg*2) ? 0.0_rt : 0.5_rt; + Real const owz = 1.0_rt - wz; Eya(j,k,l) = owx * owz * Eyc(jg ,kg ,lg ) + wx * owz * Eyc(jg+1,kg ,lg ) + owx * wz * Eyc(jg ,kg ,lg+1) + wx * wz * Eyc(jg+1,kg ,lg+1) + Eyf(j,k,l); #else - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real const wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real const owy = 1.0_rt - wy; Eya(j,k,l) = owx * owy * Eyc(jg ,kg ,lg) + wx * owy * Eyc(jg+1,kg ,lg) + owx * wy * Eyc(jg ,kg+1,lg) @@ -132,22 +132,22 @@ void warpx_interp_efield_y (int j, int k, int l, AMREX_GPU_DEVICE AMREX_FORCE_INLINE void warpx_interp_efield_z (int j, int k, int l, - amrex::Array4<amrex::Real> const& Eza, + amrex::Array4<amrex::Real > const& Eza, amrex::Array4<amrex::Real const> const& Ezf, amrex::Array4<amrex::Real const> const& Ezc) { using namespace amrex; - int lg = amrex::coarsen(l,2); - int kg = amrex::coarsen(k,2); - int jg = amrex::coarsen(j,2); + int const lg = amrex::coarsen(l,2); + int const kg = amrex::coarsen(k,2); + int const jg = amrex::coarsen(j,2); - Real wx = (j == jg*2) ? 0.0 : 0.5; - Real owx = 1.0-wx; + Real const wx = (j == jg*2) ? 0.0_rt : 0.5_rt; + Real const owx = 1.0_rt - wx; #if (AMREX_SPACEDIM == 3) - Real wy = (k == kg*2) ? 0.0 : 0.5; - Real owy = 1.0-wy; + Real wy = (k == kg*2) ? 0.0_rt : 0.5_rt; + Real owy = 1.0_rt - wy; Eza(j,k,l) = owx * owy * Ezc(jg ,kg ,lg ) + wx * owy * Ezc(jg+1,kg ,lg ) + owx * wy * Ezc(jg ,kg+1,lg ) @@ -158,4 +158,489 @@ void warpx_interp_efield_z (int j, int k, int l, #endif } +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bxa, + amrex::Array4<amrex::Real const> const& Bxf, + amrex::Array4<amrex::Real const> const& Bxc, + amrex::Array4<amrex::Real const> const& Bxg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * Bxg(jg ,kg ,0) + + owx * wy * Bxg(jg ,kg+1,0) + + wx * owy * Bxg(jg+1,kg ,0) + + wx * wy * Bxg(jg+1,kg+1,0); + + // interp from coarse staggered to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + Real bc = owx * owy * Bxc(jg ,kg ,0) + + owx * wy * Bxc(jg ,kg-1,0) + + wx * owy * Bxc(jg+1,kg ,0) + + wx * wy * Bxc(jg+1,kg-1,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5*(Bxf(j,k-1,0) + Bxf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * owz * Bxg(jg ,kg ,lg ) + + wx * owy * owz * Bxg(jg+1,kg ,lg ) + + owx * wy * owz * Bxg(jg ,kg+1,lg ) + + wx * wy * owz * Bxg(jg+1,kg+1,lg ) + + owx * owy * wz * Bxg(jg ,kg ,lg+1) + + wx * owy * wz * Bxg(jg+1,kg ,lg+1) + + owx * wy * wz * Bxg(jg ,kg+1,lg+1) + + wx * wy * wz * Bxg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + wz = 0.5-wz; owz = 1.0-wz; + Real bc = owx * owy * owz * Bxc(jg ,kg ,lg ) + + wx * owy * owz * Bxc(jg+1,kg ,lg ) + + owx * wy * owz * Bxc(jg ,kg-1,lg ) + + wx * wy * owz * Bxc(jg+1,kg-1,lg ) + + owx * owy * wz * Bxc(jg ,kg ,lg-1) + + wx * owy * wz * Bxc(jg+1,kg ,lg-1) + + owx * wy * wz * Bxc(jg ,kg-1,lg-1) + + wx * wy * wz * Bxc(jg+1,kg-1,lg-1); + + // interp from fine stagged to fine nodal + Real bf = 0.25*(Bxf(j,k-1,l-1) + Bxf(j,k,l-1) + Bxf(j,k-1,l) + Bxf(j,k,l)); +#endif + + Bxa(j,k,l) = bg + (bf-bc); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bya, + amrex::Array4<amrex::Real const> const& Byf, + amrex::Array4<amrex::Real const> const& Byc, + amrex::Array4<amrex::Real const> const& Byg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * Byg(jg ,kg ,0) + + owx * wy * Byg(jg ,kg+1,0) + + wx * owy * Byg(jg+1,kg ,0) + + wx * wy * Byg(jg+1,kg+1,0); + + // interp from coarse stagged (cell-centered for By) to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + wy = 0.5-wy; owy = 1.0-wy; + Real bc = owx * owy * Byc(jg ,kg ,0) + + owx * wy * Byc(jg ,kg-1,0) + + wx * owy * Byc(jg-1,kg ,0) + + wx * wy * Byc(jg-1,kg-1,0); + + // interp form fine stagger (cell-centered for By) to fine nodal + Real bf = 0.25*(Byf(j,k,0) + Byf(j-1,k,0) + Byf(j,k-1,0) + Byf(j-1,k-1,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * owz * Byg(jg ,kg ,lg ) + + wx * owy * owz * Byg(jg+1,kg ,lg ) + + owx * wy * owz * Byg(jg ,kg+1,lg ) + + wx * wy * owz * Byg(jg+1,kg+1,lg ) + + owx * owy * wz * Byg(jg ,kg ,lg+1) + + wx * owy * wz * Byg(jg+1,kg ,lg+1) + + owx * wy * wz * Byg(jg ,kg+1,lg+1) + + wx * wy * wz * Byg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + wz = 0.5-wz; owz = 1.0-wz; + Real bc = owx * owy * owz * Byc(jg ,kg ,lg ) + + wx * owy * owz * Byc(jg-1,kg ,lg ) + + owx * wy * owz * Byc(jg ,kg+1,lg ) + + wx * wy * owz * Byc(jg-1,kg+1,lg ) + + owx * owy * wz * Byc(jg ,kg ,lg-1) + + wx * owy * wz * Byc(jg-1,kg ,lg-1) + + owx * wy * wz * Byc(jg ,kg+1,lg-1) + + wx * wy * wz * Byc(jg-1,kg+1,lg-1); + + // interp from fine stagged to fine nodal + Real bf = 0.25*(Byf(j-1,k,l-1) + Byf(j,k,l-1) + Byf(j-1,k,l) + Byf(j,k,l)); + +#endif + + Bya(j,k,l) = bg + (bf-bc); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bza, + amrex::Array4<amrex::Real const> const& Bzf, + amrex::Array4<amrex::Real const> const& Bzc, + amrex::Array4<amrex::Real const> const& Bzg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * Bzg(jg ,kg ,0) + + owx * wy * Bzg(jg ,kg+1,0) + + wx * owy * Bzg(jg+1,kg ,0) + + wx * wy * Bzg(jg+1,kg+1,0); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + Real bc = owx * owy * Bzc(jg ,kg ,0) + + owx * wy * Bzc(jg ,kg+1,0) + + wx * owy * Bzc(jg-1,kg ,0) + + wx * wy * Bzc(jg-1,kg+1,0); + + // interp from fine staggered to fine nodal + Real bf = 0.5*(Bzf(j-1,k,0) + Bzf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real bg = owx * owy * owz * Bzg(jg ,kg ,lg ) + + wx * owy * owz * Bzg(jg+1,kg ,lg ) + + owx * wy * owz * Bzg(jg ,kg+1,lg ) + + wx * wy * owz * Bzg(jg+1,kg+1,lg ) + + owx * owy * wz * Bzg(jg ,kg ,lg+1) + + wx * owy * wz * Bzg(jg+1,kg ,lg+1) + + owx * wy * wz * Bzg(jg ,kg+1,lg+1) + + wx * wy * wz * Bzg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + wy = 0.5-wy; owy = 1.0-wy; + Real bc = owx * owy * owz * Bzc(jg ,kg ,lg ) + + wx * owy * owz * Bzc(jg-1,kg ,lg ) + + owx * wy * owz * Bzc(jg ,kg-1,lg ) + + wx * wy * owz * Bzc(jg-1,kg-1,lg ) + + owx * owy * wz * Bzc(jg ,kg ,lg+1) + + wx * owy * wz * Bzc(jg-1,kg ,lg+1) + + owx * wy * wz * Bzc(jg ,kg-1,lg+1) + + wx * wy * wz * Bzc(jg-1,kg-1,lg+1); + + // interp from fine stagged to fine nodal + Real bf = 0.25*(Bzf(j-1,k-1,l) + Bzf(j,k-1,l) + Bzf(j-1,k,l) + Bzf(j,k,l)); + +#endif + + Bza(j,k,l) = bg + (bf-bc); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bxa, + amrex::Array4<amrex::Real const> const& Bxf) +{ +#if (AMREX_SPACEDIM == 2) + Bxa(j,k,0) = 0.5*(Bxf(j,k-1,0) + Bxf(j,k,0)); +#else + Bxa(j,k,l) = 0.25*(Bxf(j,k-1,l-1) + Bxf(j,k,l-1) + Bxf(j,k-1,l) + Bxf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bya, + amrex::Array4<amrex::Real const> const& Byf) +{ +#if (AMREX_SPACEDIM == 2) + Bya(j,k,0) = 0.25*(Byf(j,k,0) + Byf(j-1,k,0) + Byf(j,k-1,0) + Byf(j-1,k-1,0)); +#else + Bya(j,k,l) = 0.25*(Byf(j-1,k,l-1) + Byf(j,k,l-1) + Byf(j-1,k,l) + Byf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_bfield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Bza, + amrex::Array4<amrex::Real const> const& Bzf) +{ +#if (AMREX_SPACEDIM == 2) + Bza(j,k,0) = 0.5*(Bzf(j-1,k,0) + Bzf(j,k,0)); +#else + Bza(j,k,l) = 0.25*(Bzf(j-1,k-1,l) + Bzf(j,k-1,l) + Bzf(j-1,k,l) + Bzf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Exa, + amrex::Array4<amrex::Real const> const& Exf, + amrex::Array4<amrex::Real const> const& Exc, + amrex::Array4<amrex::Real const> const& Exg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * Exg(jg ,kg ,0) + + owx * wy * Exg(jg ,kg+1,0) + + wx * owy * Exg(jg+1,kg ,0) + + wx * wy * Exg(jg+1,kg+1,0); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + Real ec = owx * owy * Exc(jg ,kg ,0) + + owx * wy * Exc(jg ,kg+1,0) + + wx * owy * Exc(jg-1,kg ,0) + + wx * wy * Exc(jg-1,kg+1,0); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Exf(j-1,k,0) + Exf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * owz * Exg(jg ,kg ,lg ) + + wx * owy * owz * Exg(jg+1,kg ,lg ) + + owx * wy * owz * Exg(jg ,kg+1,lg ) + + wx * wy * owz * Exg(jg+1,kg+1,lg ) + + owx * owy * wz * Exg(jg ,kg ,lg+1) + + wx * owy * wz * Exg(jg+1,kg ,lg+1) + + owx * wy * wz * Exg(jg ,kg+1,lg+1) + + wx * wy * wz * Exg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wx = 0.5-wx; owx = 1.0-wx; + Real ec = owx * owy * owz * Exc(jg ,kg ,lg ) + + wx * owy * owz * Exc(jg-1,kg ,lg ) + + owx * wy * owz * Exc(jg ,kg+1,lg ) + + wx * wy * owz * Exc(jg-1,kg+1,lg ) + + owx * owy * wz * Exc(jg ,kg ,lg+1) + + wx * owy * wz * Exc(jg-1,kg ,lg+1) + + owx * wy * wz * Exc(jg ,kg+1,lg+1) + + wx * wy * wz * Exc(jg-1,kg+1,lg+1); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Exf(j-1,k,l) + Exf(j,k,l)); + +#endif + + Exa(j,k,l) = eg + (ef-ec); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eya, + amrex::Array4<amrex::Real const> const& Eyf, + amrex::Array4<amrex::Real const> const& Eyc, + amrex::Array4<amrex::Real const> const& Eyg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal and coarse staggered to fine nodal + Real eg = owx * owy * (Eyg(jg ,kg ,0) + Eyc(jg ,kg ,0)) + + owx * wy * (Eyg(jg ,kg+1,0) + Eyc(jg ,kg+1,0)) + + wx * owy * (Eyg(jg+1,kg ,0) + Eyc(jg+1,kg ,0)) + + wx * wy * (Eyg(jg+1,kg+1,0) + Eyc(jg+1,kg+1,0)); + Real ec = 0.0; + + // interp from fine staggered to fine nodal + Real ef = Eyf(j,k,0); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * owz * Eyg(jg ,kg ,lg ) + + wx * owy * owz * Eyg(jg+1,kg ,lg ) + + owx * wy * owz * Eyg(jg ,kg+1,lg ) + + wx * wy * owz * Eyg(jg+1,kg+1,lg ) + + owx * owy * wz * Eyg(jg ,kg ,lg+1) + + wx * owy * wz * Eyg(jg+1,kg ,lg+1) + + owx * wy * wz * Eyg(jg ,kg+1,lg+1) + + wx * wy * wz * Eyg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + Real ec = owx * owy * owz * Eyc(jg ,kg ,lg ) + + wx * owy * owz * Eyc(jg+1,kg ,lg ) + + owx * wy * owz * Eyc(jg ,kg-1,lg ) + + wx * wy * owz * Eyc(jg+1,kg-1,lg ) + + owx * owy * wz * Eyc(jg ,kg ,lg+1) + + wx * owy * wz * Eyc(jg+1,kg ,lg+1) + + owx * wy * wz * Eyc(jg ,kg-1,lg+1) + + wx * wy * wz * Eyc(jg+1,kg-1,lg+1); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Eyf(j,k-1,l) + Eyf(j,k,l)); + +#endif + + Eya(j,k,l) = eg + (ef-ec); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eza, + amrex::Array4<amrex::Real const> const& Ezf, + amrex::Array4<amrex::Real const> const& Ezc, + amrex::Array4<amrex::Real const> const& Ezg) +{ + using namespace amrex; + + int jg = amrex::coarsen(j,2); + Real wx = (j == jg*2) ? 0.0 : 0.5; + Real owx = 1.0-wx; + + int kg = amrex::coarsen(k,2); + Real wy = (k == kg*2) ? 0.0 : 0.5; + Real owy = 1.0-wy; + +#if (AMREX_SPACEDIM == 2) + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * Ezg(jg ,kg ,0) + + owx * wy * Ezg(jg ,kg+1,0) + + wx * owy * Ezg(jg+1,kg ,0) + + wx * wy * Ezg(jg+1,kg+1,0); + + // interp from coarse stagged to fine nodal + wy = 0.5-wy; owy = 1.0-wy; + Real ec = owx * owy * Ezc(jg ,kg ,0) + + owx * wy * Ezc(jg ,kg-1,0) + + wx * owy * Ezc(jg+1,kg ,0) + + wx * wy * Ezc(jg+1,kg-1,0); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Ezf(j,k-1,0) + Ezf(j,k,0)); + +#else + + int lg = amrex::coarsen(l,2); + Real wz = (l == lg*2) ? 0.0 : 0.5; + Real owz = 1.0-wz; + + // interp from coarse nodal to fine nodal + Real eg = owx * owy * owz * Ezg(jg ,kg ,lg ) + + wx * owy * owz * Ezg(jg+1,kg ,lg ) + + owx * wy * owz * Ezg(jg ,kg+1,lg ) + + wx * wy * owz * Ezg(jg+1,kg+1,lg ) + + owx * owy * wz * Ezg(jg ,kg ,lg+1) + + wx * owy * wz * Ezg(jg+1,kg ,lg+1) + + owx * wy * wz * Ezg(jg ,kg+1,lg+1) + + wx * wy * wz * Ezg(jg+1,kg+1,lg+1); + + // interp from coarse staggered to fine nodal + wz = 0.5-wz; owz = 1.0-wz; + Real ec = owx * owy * owz * Ezc(jg ,kg ,lg ) + + wx * owy * owz * Ezc(jg+1,kg ,lg ) + + owx * wy * owz * Ezc(jg ,kg+1,lg ) + + wx * wy * owz * Ezc(jg+1,kg+1,lg ) + + owx * owy * wz * Ezc(jg ,kg ,lg-1) + + wx * owy * wz * Ezc(jg+1,kg ,lg-1) + + owx * wy * wz * Ezc(jg ,kg+1,lg-1) + + wx * wy * wz * Ezc(jg+1,kg+1,lg-1); + + // interp from fine staggered to fine nodal + Real ef = 0.5*(Ezf(j,k,l-1) + Ezf(j,k,l)); + +#endif + + Eza(j,k,l) = eg + (ef-ec); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_x (int j, int k, int l, + amrex::Array4<amrex::Real> const& Exa, + amrex::Array4<amrex::Real const> const& Exf) +{ + Exa(j,k,l) = 0.5*(Exf(j-1,k,l) + Exf(j,k,l)); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_y (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eya, + amrex::Array4<amrex::Real const> const& Eyf) +{ +#if (AMREX_SPACEDIM == 2) + Eya(j,k,0) = Eyf(j,k,0); +#else + Eya(j,k,l) = 0.5*(Eyf(j,k-1,l) + Eyf(j,k,l)); +#endif +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void warpx_interp_nd_efield_z (int j, int k, int l, + amrex::Array4<amrex::Real> const& Eza, + amrex::Array4<amrex::Real const> const& Ezf) +{ +#if (AMREX_SPACEDIM == 2) + Eza(j,k,0) = 0.5*(Ezf(j,k-1,0) + Ezf(j,k,0)); +#else + Eza(j,k,l) = 0.5*(Ezf(j,k,l-1) + Ezf(j,k,l)); +#endif +} + #endif diff --git a/Source/Parallelization/WarpXRegrid.cpp b/Source/Parallelization/WarpXRegrid.cpp index 5441755f5..2ae167283 100644 --- a/Source/Parallelization/WarpXRegrid.cpp +++ b/Source/Parallelization/WarpXRegrid.cpp @@ -91,7 +91,7 @@ WarpX::RemakeLevel (int lev, Real time, const BoxArray& ba, const DistributionMa // Aux patch - if (lev == 0) + if (lev == 0 && Bfield_aux[0][0]->ixType() == Bfield_fp[0][0]->ixType()) { for (int idim = 0; idim < 3; ++idim) { Bfield_aux[lev][idim].reset(new MultiFab(*Bfield_fp[lev][idim], amrex::make_alias, 0, Bfield_aux[lev][idim]->nComp())); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 6da0f1155..2737eb008 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -4,6 +4,9 @@ #include "ShapeFactors.H" #include <WarpX_Complex.H> +#include <AMReX_Array4.H> +#include <AMReX_REAL.H> + /* \brief Current Deposition for thread thread_num * /param xp, yp, zp : Pointer to arrays of particle positions. * \param wp : Pointer to array of particle weights. @@ -208,69 +211,71 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, const amrex::Real q, const long n_rz_azimuthal_modes) { + using namespace amrex; + // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) - const bool do_ionization = ion_lev; - const amrex::Real dxi = 1.0/dx[0]; - const amrex::Real dtsdx0 = dt*dxi; - const amrex::Real xmin = xyzmin[0]; + bool const do_ionization = ion_lev; + Real const dxi = 1.0_rt / dx[0]; + Real const dtsdx0 = dt*dxi; + Real const xmin = xyzmin[0]; #if (defined WARPX_DIM_3D) - const amrex::Real dyi = 1.0/dx[1]; - const amrex::Real dtsdy0 = dt*dyi; - const amrex::Real ymin = xyzmin[1]; + Real const dyi = 1.0_rt / dx[1]; + Real const dtsdy0 = dt*dyi; + Real const ymin = xyzmin[1]; #endif - const amrex::Real dzi = 1.0/dx[2]; - const amrex::Real dtsdz0 = dt*dzi; - const amrex::Real zmin = xyzmin[2]; + Real const dzi = 1.0_rt / dx[2]; + Real const dtsdz0 = dt*dzi; + Real const zmin = xyzmin[2]; #if (defined WARPX_DIM_3D) - const amrex::Real invdtdx = 1.0/(dt*dx[1]*dx[2]); - const amrex::Real invdtdy = 1.0/(dt*dx[0]*dx[2]); - const amrex::Real invdtdz = 1.0/(dt*dx[0]*dx[1]); + Real const invdtdx = 1.0_rt / (dt*dx[1]*dx[2]); + Real const invdtdy = 1.0_rt / (dt*dx[0]*dx[2]); + Real const invdtdz = 1.0_rt / (dt*dx[0]*dx[1]); #elif (defined WARPX_DIM_XZ) || (defined WARPX_DIM_RZ) - const amrex::Real invdtdx = 1.0/(dt*dx[2]); - const amrex::Real invdtdz = 1.0/(dt*dx[0]); - const amrex::Real invvol = 1.0/(dx[0]*dx[2]); + Real const invdtdx = 1.0_rt / (dt*dx[2]); + Real const invdtdz = 1.0_rt / (dt*dx[0]); + Real const invvol = 1.0_rt / (dx[0]*dx[2]); #endif #if (defined WARPX_DIM_RZ) - const Complex I = Complex{0., 1.}; + Complex const I = Complex{0., 1.}; #endif - const amrex::Real clightsq = 1.0/PhysConst::c/PhysConst::c; + Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); // Loop over particles and deposit into Jx_arr, Jy_arr and Jz_arr amrex::ParallelFor( np_to_depose, - [=] AMREX_GPU_DEVICE (long ip) { + [=] AMREX_GPU_DEVICE (long const ip) { // --- Get particle quantities - const amrex::Real gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + Real const gaminv = 1.0/std::sqrt(1.0 + uxp[ip]*uxp[ip]*clightsq + uyp[ip]*uyp[ip]*clightsq + uzp[ip]*uzp[ip]*clightsq); // wqx, wqy wqz are particle current in each direction - amrex::Real wq = q*wp[ip]; + Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; } - const amrex::Real wqx = wq*invdtdx; + Real const wqx = wq*invdtdx; #if (defined WARPX_DIM_3D) - const amrex::Real wqy = wq*invdtdy; + Real const wqy = wq*invdtdy; #endif - const amrex::Real wqz = wq*invdtdz; + Real const wqz = wq*invdtdz; // computes current and old position in grid units #if (defined WARPX_DIM_RZ) - const amrex::Real xp_mid = xp[ip] - 0.5*dt*uxp[ip]*gaminv; - const amrex::Real yp_mid = yp[ip] - 0.5*dt*uyp[ip]*gaminv; - const amrex::Real xp_old = xp[ip] - dt*uxp[ip]*gaminv; - const amrex::Real yp_old = yp[ip] - dt*uyp[ip]*gaminv; - const amrex::Real rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); - const amrex::Real rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid); - const amrex::Real rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); - amrex::Real costheta_new, sintheta_new; - if (rp_new > 0.) { + Real const xp_mid = xp[ip] - 0.5_rt * dt*uxp[ip]*gaminv; + Real const yp_mid = yp[ip] - 0.5_rt * dt*uyp[ip]*gaminv; + Real const xp_old = xp[ip] - dt*uxp[ip]*gaminv; + Real const yp_old = yp[ip] - dt*uyp[ip]*gaminv; + Real const rp_new = std::sqrt(xp[ip]*xp[ip] + yp[ip]*yp[ip]); + Real const rp_mid = std::sqrt(xp_mid*xp_mid + yp_mid*yp_mid); + Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); + Real costheta_new, sintheta_new; + if (rp_new > 0._rt) { costheta_new = xp[ip]/rp_new; sintheta_new = yp[ip]/rp_new; } else { @@ -278,7 +283,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, sintheta_new = 0.; } amrex::Real costheta_mid, sintheta_mid; - if (rp_mid > 0.) { + if (rp_mid > 0._rt) { costheta_mid = xp_mid/rp_mid; sintheta_mid = yp_mid/rp_mid; } else { @@ -286,7 +291,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, sintheta_mid = 0.; } amrex::Real costheta_old, sintheta_old; - if (rp_old > 0.) { + if (rp_old > 0._rt) { costheta_old = xp_old/rp_old; sintheta_old = yp_old/rp_old; } else { @@ -296,37 +301,37 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, const Complex xy_new0 = Complex{costheta_new, sintheta_new}; const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid}; const Complex xy_old0 = Complex{costheta_old, sintheta_old}; - const amrex::Real x_new = (rp_new - xmin)*dxi; - const amrex::Real x_old = (rp_old - xmin)*dxi; + Real const x_new = (rp_new - xmin)*dxi; + Real const x_old = (rp_old - xmin)*dxi; #else - const amrex::Real x_new = (xp[ip] - xmin)*dxi; - const amrex::Real x_old = x_new - dtsdx0*uxp[ip]*gaminv; + Real const x_new = (xp[ip] - xmin)*dxi; + Real const x_old = x_new - dtsdx0*uxp[ip]*gaminv; #endif #if (defined WARPX_DIM_3D) - const amrex::Real y_new = (yp[ip] - ymin)*dyi; - const amrex::Real y_old = y_new - dtsdy0*uyp[ip]*gaminv; + Real const y_new = (yp[ip] - ymin)*dyi; + Real const y_old = y_new - dtsdy0*uyp[ip]*gaminv; #endif - const amrex::Real z_new = (zp[ip] - zmin)*dzi; - const amrex::Real z_old = z_new - dtsdz0*uzp[ip]*gaminv; + Real const z_new = (zp[ip] - zmin)*dzi; + Real const z_old = z_new - dtsdz0*uzp[ip]*gaminv; #if (defined WARPX_DIM_RZ) - const amrex::Real vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; + Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; #elif (defined WARPX_DIM_XZ) - const amrex::Real vy = uyp[ip]*gaminv; + Real const vy = uyp[ip]*gaminv; #endif // Shape factor arrays // Note that there are extra values above and below // to possibly hold the factor for the old particle // which can be at a different grid location. - amrex::Real sx_new[depos_order + 3] = {0.}; - amrex::Real sx_old[depos_order + 3] = {0.}; + Real sx_new[depos_order + 3] = {0.}; + Real sx_old[depos_order + 3] = {0.}; #if (defined WARPX_DIM_3D) - amrex::Real sy_new[depos_order + 3] = {0.}; - amrex::Real sy_old[depos_order + 3] = {0.}; + Real sy_new[depos_order + 3] = {0.}; + Real sy_old[depos_order + 3] = {0.}; #endif - amrex::Real sz_new[depos_order + 3] = {0.}; - amrex::Real sz_old[depos_order + 3] = {0.}; + Real sz_new[depos_order + 3] = {0.}; + Real sz_old[depos_order + 3] = {0.}; // --- Compute shape factors // Compute shape factors for position as they are now and at old positions @@ -397,7 +402,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djr_cmplx = amrex::Real(2.)*sdxi*xy_mid; + const Complex djr_cmplx = 2._rt *sdxi*xy_mid; amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djr_cmplx.real()); amrex::Gpu::Atomic::Add( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djr_cmplx.imag()); xy_mid = xy_mid*xy_mid0; @@ -407,8 +412,8 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, } for (int k=dkl; k<=depos_order+2-dku; k++) { for (int i=dil; i<=depos_order+2-diu; i++) { - const amrex::Real sdyj = wq*vy*invvol*((sz_new[k] + 0.5*(sz_old[k] - sz_new[k]))*sx_new[i] + - (0.5*sz_new[k] + 1./3.*(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i])); + Real const sdyj = wq*vy*invvol*((sz_new[k] + 0.5_rt * (sz_old[k] - sz_new[k]))*sx_new[i] + + (0.5_rt * sz_new[k] + 1._rt / 3._rt *(sz_old[k] - sz_new[k]))*(sx_old[i] - sx_new[i])); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj); #if (defined WARPX_DIM_RZ) Complex xy_new = xy_new0; @@ -418,7 +423,7 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -amrex::Real(2.)*I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode* + const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode* (sx_new[i]*sz_new[k]*(xy_new - xy_mid) + sx_old[i]*sz_old[k]*(xy_mid - xy_old)); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); amrex::Gpu::Atomic::Add( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djt_cmplx.imag()); @@ -430,15 +435,15 @@ void doEsirkepovDepositionShapeN (const amrex::ParticleReal * const xp, } } for (int i=dil; i<=depos_order+2-diu; i++) { - amrex::Real sdzk = 0.; + Real sdzk = 0.; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5*(sx_old[i] - sx_new[i])); + sdzk += wqz*(sz_old[k] - sz_new[k])*(sx_new[i] + 0.5_rt * (sx_old[i] - sx_new[i])); amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdzk); #if (defined WARPX_DIM_RZ) Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes - const Complex djz_cmplx = amrex::Real(2.)*sdzk*xy_mid; + const Complex djz_cmplx = 2._rt * sdzk * xy_mid; amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djz_cmplx.real()); amrex::Gpu::Atomic::Add( &Jz_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode), djz_cmplx.imag()); xy_mid = xy_mid*xy_mid0; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 58546a106..30f7354d0 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -162,9 +162,9 @@ 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 nSpeciesBackTransformedDiagnostics() const {return nspecies_back_transformed_diagnostics;} + int mapSpeciesBackTransformedDiagnostics(int i) const {return map_species_back_transformed_diagnostics[i];} + int doBackTransformedDiagnostics() const {return do_back_transformed_diagnostics;} int nSpeciesDepositOnMainGrid () const { bool const onMainGrid = true; @@ -215,8 +215,8 @@ protected: #ifdef WARPX_QED // The QED engines - BreitWheelerEngine bw_engine; - QuantumSynchrotronEngine qs_engine; + std::shared_ptr<BreitWheelerEngine> shr_p_bw_engine; + std::shared_ptr<QuantumSynchrotronEngine> shr_p_qs_engine; //_______________________________ //Initialize QED engines and provides smart pointers @@ -236,12 +236,12 @@ private: void mapSpeciesProduct (); int getSpeciesID (std::string product_str); - // Number of species dumped in BoostedFrameDiagnostics - int nspecies_boosted_frame_diags = 0; - // map_species_boosted_frame_diags[i] is the species ID in - // MultiParticleContainer for 0<i<nspecies_boosted_frame_diags - std::vector<int> map_species_boosted_frame_diags; - int do_boosted_frame_diags = 0; + // Number of species dumped in BackTransformedDiagnostics + int nspecies_back_transformed_diagnostics = 0; + // map_species_back_transformed_diagnostics[i] is the species ID in + // MultiParticleContainer for 0<i<nspecies_back_transformed_diagnostics + std::vector<int> map_species_back_transformed_diagnostics; + int do_back_transformed_diagnostics = 0; // runtime parameters int nlasers = 0; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index c860d21f5..63aa500e9 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -38,14 +38,14 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) // 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; + map_species_back_transformed_diagnostics.resize(nspecies); + nspecies_back_transformed_diagnostics = 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 (pc->do_back_transformed_diagnostics){ + map_species_back_transformed_diagnostics[nspecies_back_transformed_diagnostics] = i; + do_back_transformed_diagnostics = 1; + nspecies_back_transformed_diagnostics += 1; } } } @@ -387,8 +387,8 @@ MultiParticleContainer BL_PROFILE("MultiParticleContainer::GetLabFrameData"); // Loop over particle species - for (int i = 0; i < nspecies_boosted_frame_diags; ++i){ - int isp = map_species_boosted_frame_diags[i]; + for (int i = 0; i < nspecies_back_transformed_diagnostics; ++i){ + int isp = map_species_back_transformed_diagnostics[i]; WarpXParticleContainer* pc = allcontainers[isp].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); @@ -603,9 +603,9 @@ namespace } // --- product runtime attribs GpuArray<ParticleReal*,6> runtime_attribs_product; - bool do_boosted_product = WarpX::do_boosted_frame_diagnostic - && pc_product->DoBoostedFrameDiags(); - if (do_boosted_product) { + bool do_back_transformed_product = WarpX::do_back_transformed_diagnostics + && pc_product->doBackTransformedDiagnostics(); + if (do_back_transformed_product) { std::map<std::string, int> comps_product = pc_product->getParticleComps(); runtime_attribs_product[0] = soa_product.GetRealData(comps_product[ "xold"]).data() + np_product_old; runtime_attribs_product[1] = soa_product.GetRealData(comps_product[ "yold"]).data() + np_product_old; @@ -652,7 +652,7 @@ namespace // Update xold etc. if boosted frame diagnostics required // for product species. Fill runtime attribs with a copy of // current properties (xold = x etc.). - if (do_boosted_product) { + if (do_back_transformed_product) { runtime_attribs_product[0][ip] = p_source.pos(0); runtime_attribs_product[1][ip] = p_source.pos(1); runtime_attribs_product[2][ip] = p_source.pos(2); @@ -736,14 +736,17 @@ MultiParticleContainer::doFieldIonization () #ifdef WARPX_QED void MultiParticleContainer::InitQED () { + shr_p_qs_engine = std::make_shared<QuantumSynchrotronEngine>(); + shr_p_bw_engine = std::make_shared<BreitWheelerEngine>(); + for (auto& pc : allcontainers) { if(pc->has_quantum_sync()){ pc->set_quantum_sync_engine_ptr - (std::make_shared<QuantumSynchrotronEngine>(qs_engine)); + (shr_p_qs_engine); } if(pc->has_breit_wheeler()){ pc->set_breit_wheeler_engine_ptr - (std::make_shared<BreitWheelerEngine>(bw_engine)); + (shr_p_bw_engine); } } } diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 3c70a957f..612da01ca 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -31,7 +31,7 @@ PhotonParticleContainer::PhotonParticleContainer (AmrCore* amr_core, int ispecie pp.query("do_qed_breit_wheeler", do_qed_breit_wheeler); //Check for processes which do not make sense for photons - bool test_quantum_sync; + bool test_quantum_sync = false; pp.query("do_qed_quantum_sync", test_quantum_sync); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( test_quantum_sync == 0, @@ -45,9 +45,8 @@ void PhotonParticleContainer::InitData() { AddParticles(0); // Note - add on level 0 - if (maxLevel() > 0) { - Redistribute(); // We then redistribute - } + Redistribute(); // We then redistribute + } void @@ -74,7 +73,7 @@ PhotonParticleContainer::PushPX(WarpXParIter& pti, const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { copy_attribs(pti, x, y, z); } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 15394fcff..6cfe20dd9 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -43,7 +43,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp 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("do_back_transformed_diagnostics", do_back_transformed_diagnostics); pp.query("do_field_ionization", do_field_ionization); @@ -86,7 +86,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //variable to set plot_flags size int plot_flag_size = PIdx::nattribs; - if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if(WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) plot_flag_size += 6; #ifdef WARPX_QED @@ -441,13 +441,13 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) for (int dir=0; dir<AMREX_SPACEDIM; dir++) { if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, Real(0.)) * dx[dir]); + overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); } else { no_overlap = true; break; } if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, Real(0.)) * dx[dir]); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); } else { no_overlap = true; break; } @@ -1078,7 +1078,7 @@ PhysicalParticleContainer::Evolve (int lev, bool has_buffer = cEx || cjx; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics) { for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { @@ -1427,9 +1427,9 @@ PhysicalParticleContainer::SplitParticles(int lev) // before splitting results in a uniform distribution after splitting const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); - amrex::Vector<amrex::Real> split_offset = {dx[0]/2./ppc_nd[0], - dx[1]/2./ppc_nd[1], - dx[2]/2./ppc_nd[2]}; + amrex::Vector<Real> split_offset = {dx[0]/2._rt/ppc_nd[0], + dx[1]/2._rt/ppc_nd[1], + dx[2]/2._rt/ppc_nd[2]}; // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); @@ -1585,7 +1585,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, const ParticleReal* const AMREX_RESTRICT By = attribs[PIdx::By].dataPtr(); const ParticleReal* const AMREX_RESTRICT Bz = attribs[PIdx::Bz].dataPtr(); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags && (a_dt_type!=DtType::SecondHalf)) + if (WarpX::do_back_transformed_diagnostics && do_back_transformed_diagnostics && (a_dt_type!=DtType::SecondHalf)) { copy_attribs(pti, x, y, z); } @@ -1842,7 +1842,7 @@ 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); + AMREX_ALWAYS_ASSERT(do_back_transformed_diagnostics == 1); const int nlevs = std::max(0, finestLevel()+1); diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H index 0abedd258..0bc0f5d01 100644 --- a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H +++ b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H @@ -23,8 +23,7 @@ void UpdateMomentumBorisWithRadiationReaction( const amrex::Real uy_old = uy; const amrex::Real uz_old = uz; - //Useful constants - const amrex::Real inv_dt = 1.0/dt; + //Useful constant constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c); //Call to regular Boris pusher diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 879f4782e..8930cb176 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -273,7 +273,7 @@ public: AddIntComp(comm); } - int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + int doBackTransformedDiagnostics () const { return do_back_transformed_diagnostics; } virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector<int>& ionization_mask) @@ -316,7 +316,7 @@ protected: amrex::Gpu::ManagedVector<amrex::Real> adk_exp_prefactor; std::string physical_element; - int do_boosted_frame_diags = 1; + int do_back_transformed_diagnostics = 1; #ifdef WARPX_QED bool do_qed = false; diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index 269171a5f..7d26e7af5 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -34,9 +34,9 @@ struct ChargeDepositionAlgo { }; struct GatheringAlgo { - // Only the Standard algorithm is implemented enum { - Standard = 0 + EnergyConserving = 0, + MomentumConserving }; }; diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp index be9cdd118..4b66a0809 100644 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -34,8 +34,9 @@ const std::map<std::string, int> charge_deposition_algo_to_int = { }; const std::map<std::string, int> gathering_algo_to_int = { - {"standard", GatheringAlgo::Standard }, - {"default", GatheringAlgo::Standard } + {"energy-conserving", GatheringAlgo::EnergyConserving }, + {"momentum-conserving", GatheringAlgo::MomentumConserving }, + {"default", GatheringAlgo::EnergyConserving } }; diff --git a/Source/Utils/WarpXTagging.cpp b/Source/Utils/WarpXTagging.cpp index 8ea3211a3..91bb802e8 100644 --- a/Source/Utils/WarpXTagging.cpp +++ b/Source/Utils/WarpXTagging.cpp @@ -22,9 +22,9 @@ WarpX::ErrorEst (int lev, TagBoxArray& tags, Real time, int /*ngrow*/) for (BoxIterator bi(bx); bi.ok(); ++bi) { const IntVect& cell = bi(); - RealVect pos {AMREX_D_DECL((cell[0]+0.5)*dx[0]+problo[0], - (cell[1]+0.5)*dx[1]+problo[1], - (cell[2]+0.5)*dx[2]+problo[2])}; + RealVect pos {AMREX_D_DECL((cell[0]+0.5_rt)*dx[0]+problo[0], + (cell[1]+0.5_rt)*dx[1]+problo[1], + (cell[2]+0.5_rt)*dx[2]+problo[2])}; if (pos > fine_tag_lo && pos < fine_tag_hi) { fab(cell) = TagBox::SET; } diff --git a/Source/WarpX.H b/Source/WarpX.H index 867a33bc3..f55670cfb 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -4,7 +4,7 @@ #include "WarpXDtType.H" #include <MultiParticleContainer.H> #include <PML.H> -#include <BoostedFrameDiagnostic.H> +#include <BackTransformedDiagnostic.H> #include <BilinearFilter.H> #include <NCIGodfreyFilter.H> @@ -105,12 +105,12 @@ public: static bool serialize_ics; // Back transformation diagnostic - static bool do_boosted_frame_diagnostic; + static bool do_back_transformed_diagnostics; static std::string lab_data_directory; static int num_snapshots_lab; static amrex::Real dt_snapshots_lab; - static bool do_boosted_frame_fields; - static bool do_boosted_frame_particles; + static bool do_back_transformed_fields; + static bool do_back_transformed_particles; // Boosted frame parameters static amrex::Real gamma_boost; @@ -209,6 +209,8 @@ public: // This function does aux(lev) = fp(lev) + I(aux(lev-1)-cp(lev)). // Caller must make sure fp and cp have ghost cells filled. void UpdateAuxilaryData (); + void UpdateAuxilaryDataStagToNodal (); + void UpdateAuxilaryDataSameType (); // Fill boundary cells including coarse/fine boundaries void FillBoundaryB (); @@ -481,7 +483,7 @@ private: std::unique_ptr<MultiParticleContainer> mypc; // Boosted Frame Diagnostics - std::unique_ptr<BoostedFrameDiagnostic> myBFD; + std::unique_ptr<BackTransformedDiagnostic> myBFD; // // Fields: First array for level, second for direction diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 53c3254f6..eef033236 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -64,12 +64,12 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; -bool WarpX::do_boosted_frame_diagnostic = false; +bool WarpX::do_back_transformed_diagnostics = 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; -bool WarpX::do_boosted_frame_particles = true; +bool WarpX::do_back_transformed_fields = true; +bool WarpX::do_back_transformed_particles = true; int WarpX::num_slice_snapshots_lab = 0; Real WarpX::dt_slice_snapshots_lab; @@ -171,7 +171,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } - do_boosted_frame_particles = mypc->doBoostedFrameDiags(); + do_back_transformed_particles = mypc->doBackTransformedDiagnostics(); Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); @@ -330,8 +330,8 @@ WarpX::ReadParameters () moving_window_v *= PhysConst::c; } - pp.query("do_boosted_frame_diagnostic", do_boosted_frame_diagnostic); - if (do_boosted_frame_diagnostic) { + pp.query("do_back_transformed_diagnostics", do_back_transformed_diagnostics); + if (do_back_transformed_diagnostics) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boosted frame diagnostic."); @@ -359,7 +359,7 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); - pp.query("do_boosted_frame_fields", do_boosted_frame_fields); + pp.query("do_back_transformed_fields", do_back_transformed_fields); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); @@ -547,9 +547,13 @@ WarpX::ReadParameters () 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"); + field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); + if (field_gathering_algo == GatheringAlgo::MomentumConserving) { + // Use same shape factors in all directions, for gathering + l_lower_order_in_v = false; + } } #ifdef WARPX_USE_PSATD @@ -602,7 +606,7 @@ WarpX::ReadParameters () } } - if (do_boosted_frame_diagnostic) { + if (do_back_transformed_diagnostics) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(gamma_boost > 1.0, "gamma_boost must be > 1 to use the boost frame diagnostic"); pp.query("num_slice_snapshots_lab", num_slice_snapshots_lab); @@ -817,16 +821,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm ncomps = n_rz_azimuthal_modes*2 - 1; #endif + bool aux_is_nodal = (field_gathering_algo == GatheringAlgo::MomentumConserving); + IntVect ngextra(static_cast<int>(aux_is_nodal and !do_nodal)); + // // The fine patch // - Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE)); - Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE)); - Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE)); + Bfield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Bx_nodal_flag),dm,ncomps,ngE+ngextra)); + Bfield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,By_nodal_flag),dm,ncomps,ngE+ngextra)); + Bfield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Bz_nodal_flag),dm,ncomps,ngE+ngextra)); - Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE)); - Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE)); - Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE)); + Efield_fp[lev][0].reset( new MultiFab(amrex::convert(ba,Ex_nodal_flag),dm,ncomps,ngE+ngextra)); + Efield_fp[lev][1].reset( new MultiFab(amrex::convert(ba,Ey_nodal_flag),dm,ncomps,ngE+ngextra)); + Efield_fp[lev][2].reset( new MultiFab(amrex::convert(ba,Ez_nodal_flag),dm,ncomps,ngE+ngextra)); current_fp[lev][0].reset( new MultiFab(amrex::convert(ba,jx_nodal_flag),dm,ncomps,ngJ)); current_fp[lev][1].reset( new MultiFab(amrex::convert(ba,jy_nodal_flag),dm,ncomps,ngJ)); @@ -873,7 +880,19 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm // // The Aux patch (i.e., the full solution) // - if (lev == 0) + if (aux_is_nodal and !do_nodal) + { + // Create aux multifabs on Nodal Box Array + BoxArray const nba = amrex::convert(ba,IntVect::TheNodeVector()); + Bfield_aux[lev][0].reset( new MultiFab(nba,dm,ncomps,ngE)); + Bfield_aux[lev][1].reset( new MultiFab(nba,dm,ncomps,ngE)); + Bfield_aux[lev][2].reset( new MultiFab(nba,dm,ncomps,ngE)); + + Efield_aux[lev][0].reset( new MultiFab(nba,dm,ncomps,ngE)); + Efield_aux[lev][1].reset( new MultiFab(nba,dm,ncomps,ngE)); + Efield_aux[lev][2].reset( new MultiFab(nba,dm,ncomps,ngE)); + } + else if (lev == 0) { for (int idir = 0; idir < 3; ++idir) { Efield_aux[lev][idir].reset(new MultiFab(*Efield_fp[lev][idir], amrex::make_alias, 0, ncomps)); @@ -954,15 +973,25 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm cba.coarsen(refRatio(lev-1)); if (n_field_gather_buffer > 0 || mypc->nSpeciesGatherFromMainGrid() > 0) { - // Create the MultiFabs for B - Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); - Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); - Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); - - // Create the MultiFabs for E - Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); - Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); - Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); + if (aux_is_nodal) { + BoxArray const& cnba = amrex::convert(cba,IntVect::TheNodeVector()); + Bfield_cax[lev][0].reset( new MultiFab(cnba,dm,ncomps,ngE)); + Bfield_cax[lev][1].reset( new MultiFab(cnba,dm,ncomps,ngE)); + Bfield_cax[lev][2].reset( new MultiFab(cnba,dm,ncomps,ngE)); + Efield_cax[lev][0].reset( new MultiFab(cnba,dm,ncomps,ngE)); + Efield_cax[lev][1].reset( new MultiFab(cnba,dm,ncomps,ngE)); + Efield_cax[lev][2].reset( new MultiFab(cnba,dm,ncomps,ngE)); + } else { + // Create the MultiFabs for B + Bfield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Bx_nodal_flag),dm,ncomps,ngE)); + Bfield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,By_nodal_flag),dm,ncomps,ngE)); + Bfield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Bz_nodal_flag),dm,ncomps,ngE)); + + // Create the MultiFabs for E + Efield_cax[lev][0].reset( new MultiFab(amrex::convert(cba,Ex_nodal_flag),dm,ncomps,ngE)); + Efield_cax[lev][1].reset( new MultiFab(amrex::convert(cba,Ey_nodal_flag),dm,ncomps,ngE)); + Efield_cax[lev][2].reset( new MultiFab(amrex::convert(cba,Ez_nodal_flag),dm,ncomps,ngE)); + } gather_buffer_masks[lev].reset( new iMultiFab(ba, dm, ncomps, 1) ); // Gather buffer masks have 1 ghost cell, because of the fact |