diff --git a/Examples/Tests/repelling_particles/analysis_repelling.py b/Examples/Tests/repelling_particles/analysis_repelling.py index ebff010fc40..bda3d74d274 100755 --- a/Examples/Tests/repelling_particles/analysis_repelling.py +++ b/Examples/Tests/repelling_particles/analysis_repelling.py @@ -36,7 +36,7 @@ # Check plotfile name specified in command line last_filename = sys.argv[1] -filename_radical = re.findall('(.*?)\d+/*$', last_filename)[0] +filename_radical = re.findall(r'(.*?)\d+/*$', last_filename)[0] # Loop through files, and extract the position and velocity of both particles x1 = [] @@ -48,10 +48,10 @@ ds = yt.load(filename) ad = ds.all_data() - x1.append( float(ad[('electron1','particle_position_x')]) ) - x2.append( float(ad[('electron2','particle_position_x')]) ) - beta1.append( float(ad[('electron1','particle_momentum_x')])/(m_e*c) ) - beta2.append( float(ad[('electron2','particle_momentum_x')])/(m_e*c) ) + x1.append( float(ad[('electron1','particle_position_x')][0]) ) + x2.append( float(ad[('electron2','particle_position_x')][0]) ) + beta1.append( float(ad[('electron1','particle_momentum_x')][0])/(m_e*c) ) + beta2.append( float(ad[('electron2','particle_momentum_x')][0])/(m_e*c) ) # Convert to numpy array x1 = np.array(x1) diff --git a/Regression/Checksum/benchmarks_json/BTD_rz.json b/Regression/Checksum/benchmarks_json/BTD_rz.json index 01e4c687292..cf83ce01346 100644 --- a/Regression/Checksum/benchmarks_json/BTD_rz.json +++ b/Regression/Checksum/benchmarks_json/BTD_rz.json @@ -1,13 +1,13 @@ { "lev=0": { - "Br": 1.8705552264208163e-08, - "Bt": 2380179.5677922587, - "Bz": 2.4079077116116535e-08, + "Br": 1.8705569155952808e-08, + "Bt": 2380179.567792259, + "Bz": 2.4079063431898616e-08, "Er": 497571119514841.25, - "Et": 7.048225464720808, - "Ez": 137058870936728.28, + "Et": 7.048223219552, + "Ez": 137058870936728.23, "jr": 0.0, "jt": 0.0, "jz": 0.0 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/PEC_particle.json b/Regression/Checksum/benchmarks_json/PEC_particle.json index 53253cebb0a..0c502a6d615 100644 --- a/Regression/Checksum/benchmarks_json/PEC_particle.json +++ b/Regression/Checksum/benchmarks_json/PEC_particle.json @@ -1,31 +1,31 @@ { "lev=0": { - "Bx": 6.52080713617892e-07, - "By": 5.323555640562078e-18, - "Bz": 1.1834601056649145e-06, - "Ex": 300.06725745121753, - "Ey": 320.82687281173185, - "Ez": 185.62208934593673, - "jx": 1.857782359527958e-06, - "jy": 257766.55567863808, + "Bx": 6.5208071361789e-07, + "By": 5.323579900863346e-18, + "Bz": 1.1834601056649094e-06, + "Ex": 300.0672574512165, + "Ey": 320.82687281173, + "Ez": 185.62208934593605, + "jx": 1.8577823595279576e-06, + "jy": 257766.5556786359, "jz": 0.0 }, + "electron": { + "particle_momentum_x": 1.1142179380730402e-32, + "particle_momentum_y": 4.5160903214571184e-36, + "particle_momentum_z": 5.735409609139057e-50, + "particle_position_x": 3.199800000000006e-05, + "particle_position_y": 6.267072536608953e-23, + "particle_position_z": 3.3938629027782103e-37, + "particle_weight": 1.0 + }, "proton": { - "particle_momentum_x": 1.0975605406855061e-33, + "particle_momentum_x": 1.097560540685506e-33, "particle_momentum_y": 1.002878875615427e-18, - "particle_momentum_z": 5.084424741580567e-51, + "particle_momentum_z": 8.484417024161229e-51, "particle_position_x": 3.1998e-05, "particle_position_y": 6.572670690061996e-06, - "particle_position_z": 5.681444715216788e-38, - "particle_weight": 1.0 - }, - "electron": { - "particle_momentum_x": 1.1142179380730476e-32, - "particle_momentum_y": 4.516090321457367e-36, - "particle_momentum_z": 2.2354046360787e-50, - "particle_position_x": 3.199800000000006e-05, - "particle_position_y": 6.267072536609136e-23, - "particle_position_z": 1.4145775327446577e-37, + "particle_position_z": 4.372139867361774e-39, "particle_weight": 1.0 } } diff --git a/Regression/Checksum/benchmarks_json/Python_ionization.json b/Regression/Checksum/benchmarks_json/Python_ionization.json index a5e65fcf765..7918a30244d 100644 --- a/Regression/Checksum/benchmarks_json/Python_ionization.json +++ b/Regression/Checksum/benchmarks_json/Python_ionization.json @@ -3,30 +3,30 @@ "Bx": 0.0, "By": 26296568.434868, "Bz": 0.0, - "Ex": 7878103122971888.0, + "Ex": 7878103122971890.0, "Ey": 0.0, - "Ez": 3027.995293554466, - "jx": 1.2111358330750162e+16, + "Ez": 3027.738911163427, + "jx": 1.2111358330750164e+16, "jy": 0.0, - "jz": 1.3483401471475687e-07 + "jz": 1.3523766702993914e-07 + }, + "electrons": { + "particle_momentum_x": 4.4237309006289776e-18, + "particle_momentum_y": 0.0, + "particle_momentum_z": 2.643146520891287e-18, + "particle_orig_z": 0.4306994012391905, + "particle_position_x": 0.11013574385450274, + "particle_position_y": 0.6415447480129455, + "particle_weight": 3.4456249999999997e-10 }, "ions": { "particle_ionizationLevel": 72897.0, - "particle_momentum_x": 1.76132401934254e-18, + "particle_momentum_x": 1.7613240193425407e-18, "particle_momentum_y": 0.0, - "particle_momentum_z": 3.644887053263054e-23, + "particle_momentum_z": 3.6448870533515125e-23, "particle_orig_z": 0.128, "particle_position_x": 0.03200001189420337, "particle_position_y": 0.1280000046901387, "particle_weight": 9.999999999999999e-11 - }, - "electrons": { - "particle_momentum_x": 4.4206237143449475e-18, - "particle_momentum_y": 0.0, - "particle_momentum_z": 2.6361297302081026e-18, - "particle_orig_z": 0.4305565137391907, - "particle_position_x": 0.11009154442846772, - "particle_position_y": 0.6414658436421568, - "particle_weight": 3.4450781249999996e-10 } } diff --git a/Regression/Checksum/benchmarks_json/RepellingParticles.json b/Regression/Checksum/benchmarks_json/RepellingParticles.json index 57b36cdc90c..0331363053a 100644 --- a/Regression/Checksum/benchmarks_json/RepellingParticles.json +++ b/Regression/Checksum/benchmarks_json/RepellingParticles.json @@ -1,32 +1,32 @@ { + "lev=0": { + "Bx": 0.0, + "By": 10603.681421961119, + "Bz": 0.0, + "Ex": 12084698236273.04, + "Ey": 0.0, + "Ez": 16282812598779.527, + "F": 21124.39775810062, + "divE": 1.495787831639977e+18, + "jx": 496020344341010.9, + "jy": 0.0, + "jz": 10.641859867595102, + "rho": 6408706.535999998 + }, "electron1": { - "particle_momentum_x": 7.303686586478214e-23, + "particle_momentum_x": 7.303686586478215e-23, "particle_momentum_y": 0.0, - "particle_momentum_z": 1.556862195127266e-36, - "particle_position_x": 1.293487226551532e-05, - "particle_position_y": 1.777889117487881e-19, + "particle_momentum_z": 1.552083282915395e-36, + "particle_position_x": 1.293487226551533e-05, + "particle_position_y": 1.7783750855614515e-19, "particle_weight": 5000000000000.0 }, "electron2": { - "particle_momentum_x": 7.303686586477884e-23, + "particle_momentum_x": 7.303686586477887e-23, "particle_momentum_y": 0.0, - "particle_momentum_z": 1.567350238825104e-36, + "particle_momentum_z": 1.5731022908002534e-36, "particle_position_x": 1.293487226551495e-05, - "particle_position_y": 1.797289839952283e-19, + "particle_position_y": 1.7989574553330234e-19, "particle_weight": 5000000000000.0 - }, - "lev=0": { - "Bx": 0.0, - "By": 10603.68142196111, - "Bz": 0.0, - "Ex": 12084698236273.04, - "Ey": 0.0, - "Ez": 16282812598779.53, - "F": 21124.39775810061, - "divE": 1.495787831639981e+18, - "jx": 496020344341011.4, - "jy": 0.0, - "jz": 10.93513947024536, - "rho": 6408706.535999998 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/ionization_boost.json b/Regression/Checksum/benchmarks_json/ionization_boost.json index 35b2db84a1a..dbde29d19f3 100644 --- a/Regression/Checksum/benchmarks_json/ionization_boost.json +++ b/Regression/Checksum/benchmarks_json/ionization_boost.json @@ -1,30 +1,30 @@ { "lev=0": { "Bx": 0.0, - "By": 18263123.342891, + "By": 18263123.342890993, "Bz": 0.0, - "Ex": 5472992180428804.0, + "Ex": 5472992180428803.0, "Ey": 0.0, - "Ez": 922.5589707939612, - "jx": 12440856508004.96, + "Ez": 922.2798874201462, + "jx": 12440856508004.957, "jy": 0.0, - "jz": 78616.0000011086 - }, - "electrons": { - "particle_momentum_x": 2.1386770170623736e-17, - "particle_momentum_y": 0.0, - "particle_momentum_z": 1.7287241262743654e-17, - "particle_position_x": 0.11064981928849067, - "particle_position_y": 1.7121826057017473, - "particle_weight": 3.0755014319061045e-09 + "jz": 78616.00000110897 }, "ions": { "particle_ionizationLevel": 52741.0, - "particle_momentum_x": 3.630619728387593e-18, + "particle_momentum_x": 3.630619728387586e-18, "particle_momentum_y": 0.0, "particle_momentum_z": 1.0432995297946715e-13, "particle_position_x": 0.021439968272741083, "particle_position_y": 0.4742770090248555, "particle_weight": 5.000948082142308e-10 + }, + "electrons": { + "particle_momentum_x": 2.138486151076505e-17, + "particle_momentum_y": 0.0, + "particle_momentum_z": 1.7285619821244245e-17, + "particle_position_x": 0.11066682695880013, + "particle_position_y": 1.712069774427624, + "particle_weight": 3.0755014319061037e-09 } -} \ No newline at end of file +} diff --git a/Regression/Checksum/benchmarks_json/ionization_lab.json b/Regression/Checksum/benchmarks_json/ionization_lab.json index 82984141b59..99d93020532 100644 --- a/Regression/Checksum/benchmarks_json/ionization_lab.json +++ b/Regression/Checksum/benchmarks_json/ionization_lab.json @@ -3,29 +3,29 @@ "Bx": 0.0, "By": 26296568.434868, "Bz": 0.0, - "Ex": 7878103122971888.0, + "Ex": 7878103122971890.0, "Ey": 0.0, - "Ez": 3027.995293554467, - "jx": 1.2111358330750162e+16, + "Ez": 3027.7389111634266, + "jx": 1.2111358330750164e+16, "jy": 0.0, - "jz": 1.3575651149270143e-07 + "jz": 1.355962687103163e-07 + }, + "electrons": { + "particle_momentum_x": 4.4195594812814575e-18, + "particle_momentum_y": 0.0, + "particle_momentum_z": 2.6520195422151263e-18, + "particle_orig_z": 0.43004993872230357, + "particle_position_x": 0.11017033205040265, + "particle_position_y": 0.6412498318709239, + "particle_weight": 3.442265625e-10 }, "ions": { "particle_ionizationLevel": 72897.0, - "particle_momentum_x": 1.7613240052056494e-18, + "particle_momentum_x": 1.761324005205651e-18, "particle_momentum_y": 0.0, - "particle_momentum_z": 3.6186967375178554e-23, + "particle_momentum_z": 3.618696737610014e-23, "particle_position_x": 0.0320000118071032, "particle_position_y": 0.12800000462171646, "particle_weight": 9.999999999999999e-11 - }, - "electrons": { - "particle_momentum_x": 4.428135211584547e-18, - "particle_momentum_y": 0.0, - "particle_momentum_z": 2.6627518442038486e-18, - "particle_orig_z": 0.43043207622230534, - "particle_position_x": 0.11023346490802505, - "particle_position_y": 0.6417814148540429, - "particle_weight": 3.44453125e-10 } -} \ No newline at end of file +} diff --git a/Source/AcceleratorLattice/LatticeElementFinder.cpp b/Source/AcceleratorLattice/LatticeElementFinder.cpp index ba2e51715fb..ec784049760 100644 --- a/Source/AcceleratorLattice/LatticeElementFinder.cpp +++ b/Source/AcceleratorLattice/LatticeElementFinder.cpp @@ -60,7 +60,7 @@ LatticeElementFinder::UpdateIndices (int const lev, amrex::MFIter const& a_mfi, // Note that the current box is used since the box may have been updated since // the initialization in InitElementFinder. const amrex::Box box = a_mfi.tilebox(); - m_zmin = WarpX::LowerCorner(box, lev, 0._rt)[2]; + m_zmin = WarpX::LowerCorner(box, lev, 0._rt).z; m_time = warpx.gett_new(lev); if (accelerator_lattice.h_quad.nelements > 0) { diff --git a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp index a5e0004325a..2b063b99a15 100644 --- a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp +++ b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp @@ -201,8 +201,8 @@ WarpX::ApplyFieldBoundaryOnAxis (amrex::MultiFab* Er, amrex::MultiFab* Et, amrex // Lower corner of tile box physical domain // Note that this is done before the tilebox.grow so that // these do not include the guard cells. - const std::array& xyzmin = LowerCorner(tilebox, lev, 0._rt); - const amrex::Real rmin = xyzmin[0]; + const amrex::XDim3 xyzmin = LowerCorner(tilebox, lev, 0._rt); + const amrex::Real rmin = xyzmin.x; // Skip blocks that don't touch the axis if (rmin > 0._rt) { continue; } diff --git a/Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp b/Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp index 434ab1d7949..1b0c74cf7fa 100644 --- a/Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp +++ b/Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp @@ -441,8 +441,7 @@ void ColliderRelevant::ComputeDiags (int step) const int lev = 0; // define variables in preparation for field gathering - const std::array& dx = WarpX::CellSize(std::max(lev, 0)); - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(lev, 0)); const amrex::MultiFab & Ex = warpx.getField(FieldType::Efield_aux, lev,0); const amrex::MultiFab & Ey = warpx.getField(FieldType::Efield_aux, lev,1); const amrex::MultiFab & Ez = warpx.getField(FieldType::Efield_aux, lev,2); @@ -478,8 +477,7 @@ void ColliderRelevant::ComputeDiags (int step) amrex::Box box = pti.tilebox(); box.grow(ngEB); const amrex::Dim3 lo = amrex::lbound(box); - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt); const amrex::Array4 & ex_arr = Ex[pti].array(); const amrex::Array4 & ey_arr = Ey[pti].array(); const amrex::Array4 & ez_arr = Ez[pti].array(); @@ -515,7 +513,7 @@ void ColliderRelevant::ComputeDiags (int step) ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, + dinv, xyzmin, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); // compute chi amrex::Real chi = 0.0_rt; diff --git a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp index 40ef1a088e6..9c0fac101a2 100644 --- a/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldEnergy.cpp @@ -178,10 +178,10 @@ FieldEnergy::ComputeNorm2RZ(const amrex::MultiFab& field, const int lev) amrex::Box tb = convert(tilebox, field.ixType().toIntVect()); // Lower corner of tile box physical domain - const std::array& xyzmin = WarpX::LowerCorner(tilebox, lev, 0._rt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, lev, 0._rt); const Dim3 lo = lbound(tilebox); const Dim3 hi = ubound(tilebox); - const Real rmin = xyzmin[0] + (tb.ixType().nodeCentered(0) ? 0._rt : 0.5_rt*dr); + const Real rmin = xyzmin.x + (tb.ixType().nodeCentered(0) ? 0._rt : 0.5_rt*dr); const int irmin = lo.x; const int irmax = hi.x; diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index 1fc3b957ec9..1113075f4cd 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -489,11 +489,8 @@ void FieldProbe::ComputeDiags (int step) auto * const AMREX_RESTRICT idcpu = pti.GetStructOfArrays().GetIdCPUData().data(); - const auto &xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - const std::array &dx = WarpX::CellSize(lev); - - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt); + const amrex::XDim3 dinv = WarpX::InvCellSize(lev); const Dim3 lo = lbound(box); // Temporarily defining modes and interp outside ParallelFor to avoid GPU compilation errors. @@ -514,7 +511,7 @@ void FieldProbe::ComputeDiags (int step) doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, arrEx, arrEy, arrEz, arrBx, arrBy, arrBz, Extype, Eytype, Eztype, Bxtype, Bytype, Bztype, - dx_arr, xyzmin_arr, lo, temp_modes, + dinv, xyzmin, lo, temp_modes, temp_interp_order, false); //Calculate the Poynting Vector S diff --git a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp index 9adfd3f238c..811fee7a9de 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleExtrema.cpp @@ -264,8 +264,8 @@ void ParticleExtrema::ComputeDiags (int step) for (int lev = 0; lev <= level_number; ++lev) { // define variables in preparation for field gathering - const std::array& dx = WarpX::CellSize(std::max(lev, 0)); - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(lev, 0)); + const amrex::MultiFab & Ex = warpx.getField(FieldType::Efield_aux, lev,0); const amrex::MultiFab & Ey = warpx.getField(FieldType::Efield_aux, lev,1); const amrex::MultiFab & Ez = warpx.getField(FieldType::Efield_aux, lev,2); @@ -300,8 +300,7 @@ void ParticleExtrema::ComputeDiags (int step) amrex::Box box = pti.tilebox(); box.grow(ngEB); const amrex::Dim3 lo = amrex::lbound(box); - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt); const auto& ex_arr = Ex[pti].array(); const auto& ey_arr = Ey[pti].array(); const auto& ez_arr = Ez[pti].array(); @@ -337,7 +336,7 @@ void ParticleExtrema::ComputeDiags (int step) ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, + dinv, xyzmin, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); // compute chi amrex::Real chi = 0.0_rt; diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp index b26248fb2d0..d295d61b43a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp @@ -199,8 +199,8 @@ SpectralFieldDataRZ::SpectralFieldDataRZ (const int lev, #endif // Create the Hankel transformer for each box. - std::array xmax = WarpX::UpperCorner(mfi.tilebox(), lev, 0._rt); - multi_spectral_hankel_transformer[mfi] = SpectralHankelTransformer(grid_size[0], n_rz_azimuthal_modes, xmax[0]); + amrex::XDim3 const xyzmax = WarpX::UpperCorner(mfi.tilebox(), lev, 0._rt); + multi_spectral_hankel_transformer[mfi] = SpectralHankelTransformer(grid_size[0], n_rz_azimuthal_modes, xyzmax.x); } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 6631c1b88ee..602a2666b27 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -1245,11 +1245,11 @@ WarpX::ApplyInverseVolumeScalingToCurrentDensity (MultiFab* Jx, MultiFab* Jy, Mu // Lower corner of tile box physical domain // Note that this is done before the tilebox.grow so that // these do not include the guard cells. - const std::array& xyzmin = WarpX::LowerCorner(tilebox, lev, 0._rt); - const Real rmin = xyzmin[0]; - const Real rminr = xyzmin[0] + (tbr.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); - const Real rmint = xyzmin[0] + (tbt.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); - const Real rminz = xyzmin[0] + (tbz.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, lev, 0._rt); + const Real rmin = xyzmin.x; + const Real rminr = xyzmin.x + (tbr.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); + const Real rmint = xyzmin.x + (tbt.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); + const Real rminz = xyzmin.x + (tbz.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); const Dim3 lo = lbound(tilebox); const int irmin = lo.x; @@ -1416,10 +1416,10 @@ WarpX::ApplyInverseVolumeScalingToChargeDensity (MultiFab* Rho, int lev) // Lower corner of tile box physical domain // Note that this is done before the tilebox.grow so that // these do not include the guard cells. - const std::array& xyzmin = WarpX::LowerCorner(tilebox, lev, 0._rt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, lev, 0._rt); const Dim3 lo = lbound(tilebox); - const Real rmin = xyzmin[0]; - const Real rminr = xyzmin[0] + (tb.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); + const Real rmin = xyzmin.x; + const Real rminr = xyzmin.x + (tb.type(0) == NODE ? 0._rt : 0.5_rt*dx[0]); const int irmin = lo.x; const int ishift = (rminr > rmin ? 1 : 0); diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index 2912b769c57..d3d36dbfbda 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -28,8 +28,8 @@ ion_lev is a null pointer. * \param rho_fab FArrayBox of charge density, either full array or tile. * \param np_to_deposit Number of particles for which current is deposited. - * \param dx 3D cell size - * \param xyzmin Physical lower bounds of domain. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param q species charge. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry. @@ -40,37 +40,19 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, const int* ion_lev, amrex::FArrayBox& rho_fab, long np_to_deposit, - const std::array& dx, - const std::array xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, amrex::Dim3 lo, amrex::Real q, - int n_rz_azimuthal_modes) + [[maybe_unused]] int 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 dzi = 1.0_rt/dx[2]; -#if defined(WARPX_DIM_1D_Z) - const amrex::Real invvol = dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real invvol = dxi*dzi; -#elif defined(WARPX_DIM_3D) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real dyi = 1.0_rt/dx[1]; - const amrex::Real invvol = dxi*dyi*dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D) - const amrex::Real xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif - const amrex::Real zmin = xyzmin[2]; + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; amrex::Array4 const& rho_arr = rho_fab.array(); amrex::IntVect const rho_type = rho_fab.box().type(); @@ -98,19 +80,12 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, // Get particle position in grid coordinates #if defined(WARPX_DIM_RZ) const amrex::Real rp = std::sqrt(xp*xp + yp*yp); - amrex::Real costheta; - amrex::Real sintheta; - if (rp > 0.) { - costheta = xp/rp; - sintheta = yp/rp; - } else { - costheta = 1._rt; - sintheta = 0._rt; - } + const amrex::Real costheta = (rp > 0._rt ? xp/rp : 1._rt); + const amrex::Real sintheta = (rp > 0._rt ? yp/rp : 0._rt); const Complex xy0 = Complex{costheta, sintheta}; - const amrex::Real x = (rp - xmin)*dxi; + const amrex::Real x = (rp - xyzmin.x)*dinv.x; #else - const amrex::Real x = (xp - xmin)*dxi; + const amrex::Real x = (xp - xyzmin.x)*dinv.x; #endif // Compute shape factor along x @@ -125,7 +100,7 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, #endif //defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D) #if defined(WARPX_DIM_3D) // y direction - const amrex::Real y = (yp - ymin)*dyi; + const amrex::Real y = (yp - xyzmin.y)*dinv.y; amrex::Real sy[depos_order + 1] = {0._rt}; int j = 0; if (rho_type[1] == NODE) { @@ -135,7 +110,7 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, } #endif // z direction - const amrex::Real z = (zp - zmin)*dzi; + const amrex::Real z = (zp - xyzmin.z)*dinv.z; amrex::Real sz[depos_order + 1] = {0._rt}; int k = 0; if (rho_type[WARPX_ZINDEX] == NODE) { @@ -151,8 +126,7 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, &rho_arr(lo.x+k+iz, 0, 0, 0), sz[iz]*wq); } -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::AddNoRet( @@ -182,10 +156,6 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, #endif } ); - -#ifndef WARPX_DIM_RZ - amrex::ignore_unused(n_rz_azimuthal_modes); -#endif } /* \brief Perform charge deposition on a tile using shared memory @@ -198,8 +168,8 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, * \param rho_fab FArrayBox of charge density, either full array or tile. * \param ix_type * \param np_to_deposit Number of particles for which charge is deposited. - * \param dx 3D cell size - * \param xyzmin Physical lower bounds of domain. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param q species charge. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry. @@ -216,11 +186,11 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio amrex::FArrayBox& rho_fab, const amrex::IntVect& ix_type, const long np_to_deposit, - const std::array& dx, - const std::array xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3 lo, const amrex::Real q, - const int n_rz_azimuthal_modes, + [[maybe_unused]]const int n_rz_azimuthal_modes, const amrex::DenseBins& a_bins, const amrex::Box& box, const amrex::Geometry& geom, @@ -239,26 +209,8 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio // 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 dzi = 1.0_rt/dx[2]; -#if defined(WARPX_DIM_1D_Z) - const amrex::Real invvol = dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real invvol = dxi*dzi; -#elif defined(WARPX_DIM_3D) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real dyi = 1.0_rt/dx[1]; - const amrex::Real invvol = dxi*dyi*dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D) - const amrex::Real xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif - const amrex::Real zmin = xyzmin[2]; + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; amrex::Array4 const& rho_arr = rho_fab.array(); auto rho_box = rho_fab.box(); @@ -269,7 +221,6 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio // Loop over particles and deposit into rho_fab #if defined(AMREX_USE_CUDA) || defined(AMREX_USE_HIP) - const auto dxiarr = geom.InvCellSizeArray(); const auto plo = geom.ProbLoArray(); const auto domain = geom.Domain(); @@ -316,15 +267,14 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio ParticleReal xp, yp, zp; GetPosition(permutation[bin_start], xp, yp, zp); #if defined(WARPX_DIM_3D) - IntVect iv = IntVect(int(amrex::Math::floor((xp-plo[0])*dxiarr[0])), - int(amrex::Math::floor((yp-plo[1])*dxiarr[1])), - int(amrex::Math::floor((zp-plo[2])*dxiarr[2]))); + IntVect iv = IntVect(int(amrex::Math::floor((xp-plo[0])*dinv.x)), + int(amrex::Math::floor((yp-plo[1])*dinv.y)), + int(amrex::Math::floor((zp-plo[2])*dinv.z))); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - IntVect iv = IntVect( - int(amrex::Math::floor((xp-plo[0])*dxiarr[0])), - int(amrex::Math::floor((zp-plo[1])*dxiarr[1]))); + IntVect iv = IntVect(int(amrex::Math::floor((xp-plo[0])*dinv.x)), + int(amrex::Math::floor((zp-plo[1])*dinv.z))); #elif defined(WARPX_DIM_1D_Z) - IntVect iv = IntVect(int(amrex::Math::floor((zp-plo[0])*dxiarr[0]))); + IntVect iv = IntVect(int(amrex::Math::floor((zp-plo[0])*dinv.z))); #endif iv += domain.smallEnd(); getTileIndex(iv, box, true, bin_size, buffer_box); @@ -381,9 +331,9 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio sintheta = 0._rt; } const Complex xy0 = Complex{costheta, sintheta}; - const amrex::Real x = (rp - xmin)*dxi; + const amrex::Real x = (rp - xyzmin.x)*dinv.x; #else - const amrex::Real x = (xp - xmin)*dxi; + const amrex::Real x = (xp - xyzmin.x)*dinv.x; #endif // Compute shape factor along x @@ -398,7 +348,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio #endif //defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D) #if defined(WARPX_DIM_3D) // y direction - const amrex::Real y = (yp - ymin)*dyi; + const amrex::Real y = (yp - xyzmin.y)*dinv.y; amrex::Real sy[depos_order + 1] = {0._rt}; int j = 0; if (rho_type[1] == NODE) { @@ -408,7 +358,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio } #endif // z direction - const amrex::Real z = (zp - zmin)*dzi; + const amrex::Real z = (zp - xyzmin.z)*dinv.z; amrex::Real sz[depos_order + 1] = {0._rt}; int k = 0; if (rho_type[WARPX_ZINDEX] == NODE) { @@ -424,8 +374,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio &buf(lo.x+k+iz, 0, 0, 0), sz[iz]*wq); } -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) +#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) for (int iz=0; iz<=depos_order; iz++){ for (int ix=0; ix<=depos_order; ix++){ amrex::Gpu::Atomic::AddNoRet( @@ -462,10 +411,6 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio #endif // defined(AMREX_USE_CUDA) || defined(AMREX_USE_HIP) } ); - -#ifndef WARPX_DIM_RZ - amrex::ignore_unused(n_rz_azimuthal_modes); -#endif } #endif // WARPX_CHARGEDEPOSITION_H_ diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 3dce7d28534..cb56c559bc0 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -24,6 +24,7 @@ #include #include #include +#include #include /** @@ -38,17 +39,17 @@ * current positions of the particles. When different than 0, * the particle position will be temporarily modified to match * the time of the deposition. - * \param dzi, dxi, dyi The inverse cell sizes - * \param zmin, xmin, ymin The lower bounds of the domain + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param invvol The inverse volume of a grid cell * \param lo Index lower bounds of domain. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry. */ template AMREX_GPU_HOST_DEVICE AMREX_INLINE -void doDepositionShapeNKernel(const amrex::ParticleReal xp, - const amrex::ParticleReal yp, - const amrex::ParticleReal zp, +void doDepositionShapeNKernel([[maybe_unused]] const amrex::ParticleReal xp, + [[maybe_unused]] const amrex::ParticleReal yp, + [[maybe_unused]] const amrex::ParticleReal zp, const amrex::ParticleReal wq, const amrex::ParticleReal vx, const amrex::ParticleReal vy, @@ -60,28 +61,14 @@ void doDepositionShapeNKernel(const amrex::ParticleReal xp, amrex::IntVect const& jy_type, amrex::IntVect const& jz_type, const amrex::Real relative_time, - AMREX_D_DECL(const amrex::Real dzi, - const amrex::Real dxi, - const amrex::Real dyi), - AMREX_D_DECL(const amrex::Real zmin, - const amrex::Real xmin, - const amrex::Real ymin), + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Real invvol, const amrex::Dim3 lo, - const int n_rz_azimuthal_modes) + [[maybe_unused]] const int n_rz_azimuthal_modes) { using namespace amrex::literals; -#if !defined(WARPX_DIM_RZ) - amrex::ignore_unused(n_rz_azimuthal_modes); -#endif -#if defined(WARPX_DIM_1D_Z) - amrex::ignore_unused(xp, yp); -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(yp); -#endif - constexpr int zdir = WARPX_ZINDEX; constexpr int NODE = amrex::IndexType::NODE; constexpr int CELL = amrex::IndexType::CELL; @@ -92,15 +79,8 @@ void doDepositionShapeNKernel(const amrex::ParticleReal xp, const amrex::Real xpmid = xp + relative_time*vx; const amrex::Real ypmid = yp + relative_time*vy; const amrex::Real rpmid = std::sqrt(xpmid*xpmid + ypmid*ypmid); - amrex::Real costheta; - amrex::Real sintheta; - if (rpmid > 0._rt) { - costheta = xpmid/rpmid; - sintheta = ypmid/rpmid; - } else { - costheta = 1._rt; - sintheta = 0._rt; - } + const amrex::Real costheta = (rpmid > 0._rt ? xpmid/rpmid : 1._rt); + const amrex::Real sintheta = (rpmid > 0._rt ? ypmid/rpmid : 0._rt); const Complex xy0 = Complex{costheta, sintheta}; const amrex::Real wqx = wq*invvol*(+vx*costheta + vy*sintheta); const amrex::Real wqy = wq*invvol*(-vx*sintheta + vy*costheta); @@ -115,12 +95,13 @@ void doDepositionShapeNKernel(const amrex::ParticleReal xp, #if (AMREX_SPACEDIM >= 2) // x direction // Get particle position after 1/2 push back in position -#if defined(WARPX_DIM_RZ) // Keep these double to avoid bug in single precision - const double xmid = (rpmid - xmin)*dxi; +#if defined(WARPX_DIM_RZ) + const double xmid = (rpmid - xyzmin.x)*dinv.x; #else - const double xmid = ((xp - xmin) + relative_time*vx)*dxi; + const double xmid = ((xp - xyzmin.x) + relative_time*vx)*dinv.x; #endif + // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current // sx_j[xyz] shape factor along x for the centering of each current // There are only two possible centerings, node or cell centered, so at most only two shape factor @@ -155,7 +136,7 @@ void doDepositionShapeNKernel(const amrex::ParticleReal xp, #if defined(WARPX_DIM_3D) // y direction // Keep these double to avoid bug in single precision - const double ymid = ((yp - ymin) + relative_time*vy)*dyi; + const double ymid = ((yp - xyzmin.y) + relative_time*vy)*dinv.y; double sy_node[depos_order + 1] = {0.}; double sy_cell[depos_order + 1] = {0.}; int k_node = 0; @@ -182,7 +163,8 @@ void doDepositionShapeNKernel(const amrex::ParticleReal xp, // z direction // Keep these double to avoid bug in single precision - const double zmid = ((zp - zmin) + relative_time*vz)*dzi; + constexpr int zdir = WARPX_ZINDEX; + const double zmid = ((zp - xyzmin.z) + relative_time*vz)*dinv.z; double sz_node[depos_order + 1] = {0.}; double sz_cell[depos_order + 1] = {0.}; int l_node = 0; @@ -282,7 +264,7 @@ void doDepositionShapeNKernel(const amrex::ParticleReal xp, * current positions of the particles. When different than 0, * the particle position will be temporarily modified to match * the time of the deposition. - * \param dx 3D cell size + * \param dinv 3D cell size inverse * \param xyzmin Physical lower bounds of domain. * \param lo Index lower bounds of domain. * \param q species charge. @@ -300,41 +282,19 @@ void doDepositionShapeN (const GetParticlePosition& GetPosition, amrex::FArrayBox& jz_fab, long np_to_deposit, amrex::Real relative_time, - const std::array& dx, - const std::array& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, amrex::Dim3 lo, amrex::Real q, - int n_rz_azimuthal_modes) + [[maybe_unused]]int n_rz_azimuthal_modes) { using namespace amrex::literals; -#if !defined(WARPX_DIM_RZ) - amrex::ignore_unused(n_rz_azimuthal_modes); -#endif - // 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 dzi = 1.0_rt/dx[2]; -#if defined(WARPX_DIM_1D_Z) - const amrex::Real invvol = dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real invvol = dxi*dzi; -#elif defined(WARPX_DIM_3D) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real dyi = 1.0_rt/dx[1]; - const amrex::Real invvol = dxi*dyi*dzi; -#endif -#if (AMREX_SPACEDIM >= 2) - const amrex::Real xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif - const amrex::Real zmin = xyzmin[2]; + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; const amrex::Real clightsq = 1.0_rt/PhysConst::c/PhysConst::c; @@ -367,9 +327,7 @@ void doDepositionShapeN (const GetParticlePosition& GetPosition, doDepositionShapeNKernel(xp, yp, zp, wq, vx, vy, vz, jx_arr, jy_arr, jz_arr, jx_type, jy_type, jz_type, - relative_time, - AMREX_D_DECL(dzi, dxi, dyi), - AMREX_D_DECL(zmin, xmin, ymin), + relative_time, dinv, xyzmin, invvol, lo, n_rz_azimuthal_modes); } @@ -391,7 +349,7 @@ void doDepositionShapeN (const GetParticlePosition& GetPosition, ion_lev is a null pointer. * \param jx_fab,jy_fab,jz_fab FArrayBox of current density, either full array or tile. * \param np_to_deposit Number of particles for which current is deposited. - * \param dx 3D cell size + * \param dinv 3D cell size inverse * \param xyzmin Physical lower bounds of domain. * \param lo Index lower bounds of domain. * \param q species charge. @@ -411,40 +369,19 @@ void doDepositionShapeNImplicit(const GetParticlePosition& GetPosition, amrex::FArrayBox& jy_fab, amrex::FArrayBox& jz_fab, const long np_to_deposit, - const std::array& dx, - const std::array& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3 lo, const amrex::Real q, - const int n_rz_azimuthal_modes) + [[maybe_unused]]const int n_rz_azimuthal_modes) { using namespace amrex::literals; -#if !defined(WARPX_DIM_RZ) - amrex::ignore_unused(n_rz_azimuthal_modes); -#endif // 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 dzi = 1.0_rt/dx[2]; -#if defined(WARPX_DIM_1D_Z) - const amrex::Real invvol = dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real invvol = dxi*dzi; -#elif defined(WARPX_DIM_3D) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real dyi = 1.0_rt/dx[1]; - const amrex::Real invvol = dxi*dyi*dzi; -#endif -#if (AMREX_SPACEDIM >= 2) - const amrex::Real xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif - const amrex::Real zmin = xyzmin[2]; + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; amrex::Array4 const& jx_arr = jx_fab.array(); amrex::Array4 const& jy_arr = jy_fab.array(); @@ -483,9 +420,7 @@ void doDepositionShapeNImplicit(const GetParticlePosition& GetPosition, const amrex::Real relative_time = 0._rt; doDepositionShapeNKernel(xp, yp, zp, wq, vx, vy, vz, jx_arr, jy_arr, jz_arr, jx_type, jy_type, jz_type, - relative_time, - AMREX_D_DECL(dzi, dxi, dyi), - AMREX_D_DECL(zmin, xmin, ymin), + relative_time, dinv, xyzmin, invvol, lo, n_rz_azimuthal_modes); } @@ -509,7 +444,7 @@ void doDepositionShapeNImplicit(const GetParticlePosition& GetPosition, * current positions of the particles. When different than 0, * the particle position will be temporarily modified to match * the time of the deposition. - * \param dx 3D cell size + * \param dinv 3D cell size inverse * \param xyzmin Physical lower bounds of domain. * \param lo Index lower bounds of domain. * \param q species charge. @@ -527,8 +462,8 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, amrex::FArrayBox& jz_fab, long np_to_deposit, const amrex::Real relative_time, - const std::array& dx, - const std::array& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, amrex::Dim3 lo, amrex::Real q, int n_rz_azimuthal_modes, @@ -634,7 +569,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, { const unsigned int ip = permutation[ip_orig]; depositComponent(GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_buff, jx_type, - relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes, + relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes, ip, zdir, NODE, CELL, 0); } @@ -649,7 +584,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, { const unsigned int ip = permutation[ip_orig]; depositComponent(GetPosition, wp, uxp, uyp, uzp, ion_lev, jy_buff, jy_type, - relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes, + relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes, ip, zdir, NODE, CELL, 1); } @@ -664,7 +599,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, { const unsigned int ip = permutation[ip_orig]; depositComponent(GetPosition, wp, uxp, uyp, uzp, ion_lev, jz_buff, jz_type, - relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes, + relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes, ip, zdir, NODE, CELL, 2); } @@ -675,7 +610,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, // Note, you should never reach this part of the code. This funcion cannot be called unless // using HIP/CUDA, and those things are checked prior //don't use any args - ignore_unused(GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes, a_bins, box, geom, a_tbox_max_size); + ignore_unused(GetPosition, wp, uxp, uyp, uzp, ion_lev, jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, n_rz_azimuthal_modes, a_bins, box, geom, a_tbox_max_size); WARPX_ABORT_WITH_MESSAGE("Shared memory only implemented for HIP/CUDA"); #endif } @@ -698,7 +633,7 @@ void doDepositionSharedShapeN (const GetParticlePosition& GetPosition, * current positions of the particles. When different than 0, * the particle position will be temporarily modified to match * the time of the deposition. - * \param dx 3D cell size + * \param dinv 3D cell size inverse * \param xyzmin Physical lower bounds of domain. * \param lo Index lower bounds of domain. * \param q species charge. @@ -717,53 +652,28 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, long np_to_deposit, amrex::Real dt, amrex::Real relative_time, - const std::array& dx, - std::array xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, amrex::Dim3 lo, amrex::Real q, - int n_rz_azimuthal_modes) + [[maybe_unused]]int n_rz_azimuthal_modes) { using namespace amrex; using namespace amrex::literals; -#if !defined(WARPX_DIM_RZ) - ignore_unused(n_rz_azimuthal_modes); -#endif - // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) bool const do_ionization = ion_lev; -#if !defined(WARPX_DIM_1D_Z) - Real const dxi = 1.0_rt / dx[0]; -#endif -#if !defined(WARPX_DIM_1D_Z) - Real const xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - Real const dyi = 1.0_rt / dx[1]; - Real const ymin = xyzmin[1]; +#if !defined(WARPX_DIM_3D) + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; #endif - Real const dzi = 1.0_rt / dx[2]; - Real const zmin = xyzmin[2]; -#if defined(WARPX_DIM_3D) - 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) - 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]); -#elif defined(WARPX_DIM_1D_Z) - Real const invdtdz = 1.0_rt / (dt*dx[0]); - Real const invvol = 1.0_rt / (dx[2]); -#endif + amrex::XDim3 const invdtd = amrex::XDim3{(1.0_rt/dt)*dinv.y*dinv.z, + (1.0_rt/dt)*dinv.x*dinv.z, + (1.0_rt/dt)*dinv.x*dinv.y}; -#if defined(WARPX_DIM_RZ) - Complex const I = Complex{0._rt, 1._rt}; -#endif + Real constexpr clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); - Real const clightsq = 1.0_rt / ( PhysConst::c * PhysConst::c ); #if !defined(WARPX_DIM_1D_Z) Real constexpr one_third = 1.0_rt / 3.0_rt; Real constexpr one_sixth = 1.0_rt / 6.0_rt; @@ -778,7 +688,6 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, + uyp[ip]*uyp[ip]*clightsq + uzp[ip]*uzp[ip]*clightsq); - // wqx, wqy wqz are particle current in each direction Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; @@ -787,14 +696,6 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, ParticleReal xp, yp, zp; GetPosition(ip, xp, yp, zp); -#if !defined(WARPX_DIM_1D_Z) - Real const wqx = wq*invdtdx; -#endif -#if defined(WARPX_DIM_3D) - Real const wqy = wq*invdtdy; -#endif - Real const wqz = wq*invdtdz; - // computes current and old position in grid units #if defined(WARPX_DIM_RZ) Real const xp_new = xp + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv; @@ -806,51 +707,33 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new); 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_new/rp_new; - sintheta_new = yp_new/rp_new; - } else { - costheta_new = 1._rt; - sintheta_new = 0._rt; - } - amrex::Real costheta_mid, sintheta_mid; - if (rp_mid > 0._rt) { - costheta_mid = xp_mid/rp_mid; - sintheta_mid = yp_mid/rp_mid; - } else { - costheta_mid = 1._rt; - sintheta_mid = 0._rt; - } - amrex::Real costheta_old, sintheta_old; - if (rp_old > 0._rt) { - costheta_old = xp_old/rp_old; - sintheta_old = yp_old/rp_old; - } else { - costheta_old = 1._rt; - sintheta_old = 0._rt; - } + const amrex::Real costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt); + const amrex::Real sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt); + const amrex::Real costheta_new = (rp_new > 0._rt ? xp_new/rp_new : 1._rt); + const amrex::Real sintheta_new = (rp_new > 0._rt ? yp_new/rp_new : 0._rt); + const amrex::Real costheta_old = (rp_old > 0._rt ? xp_old/rp_old : 1._rt); + const amrex::Real sintheta_old = (rp_old > 0._rt ? yp_old/rp_old : 0._rt); 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}; // Keep these double to avoid bug in single precision - double const x_new = (rp_new - xmin)*dxi; - double const x_old = (rp_old - xmin)*dxi; + double const x_new = (rp_new - xyzmin.x)*dinv.x; + double const x_old = (rp_old - xyzmin.x)*dinv.x; #else #if !defined(WARPX_DIM_1D_Z) // Keep these double to avoid bug in single precision - double const x_new = (xp - xmin + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv)*dxi; - double const x_old = x_new - dt*dxi*uxp[ip]*gaminv; + double const x_new = (xp - xyzmin.x + (relative_time + 0.5_rt*dt)*uxp[ip]*gaminv)*dinv.x; + double const x_old = x_new - dt*dinv.x*uxp[ip]*gaminv; #endif #endif #if defined(WARPX_DIM_3D) // Keep these double to avoid bug in single precision - double const y_new = (yp - ymin + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv)*dyi; - double const y_old = y_new - dt*dyi*uyp[ip]*gaminv; + double const y_new = (yp - xyzmin.y + (relative_time + 0.5_rt*dt)*uyp[ip]*gaminv)*dinv.y; + double const y_old = y_new - dt*dinv.y*uyp[ip]*gaminv; #endif // Keep these double to avoid bug in single precision - double const z_new = (zp - zmin + (relative_time + 0.5_rt*dt)*uzp[ip]*gaminv)*dzi; - double const z_old = z_new - dt*dzi*uzp[ip]*gaminv; + double const z_new = (zp - xyzmin.z + (relative_time + 0.5_rt*dt)*uzp[ip]*gaminv)*dinv.z; + double const z_old = z_new - dt*dinv.z*uzp[ip]*gaminv; #if defined(WARPX_DIM_RZ) Real const vy = (-uxp[ip]*sintheta_mid + uyp[ip]*costheta_mid)*gaminv; @@ -861,6 +744,12 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, Real const vy = uyp[ip]*gaminv; #endif + // --- Compute shape factors + // Compute shape factors for position as they are now and at old positions + // [ijk]_new: leftmost grid point that the particle touches + const Compute_shape_factor< depos_order > compute_shape_factor; + const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; + // Shape factor arrays // Note that there are extra values above and below // to possibly hold the factor for the old particle @@ -869,30 +758,17 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, #if !defined(WARPX_DIM_1D_Z) double sx_new[depos_order + 3] = {0.}; double sx_old[depos_order + 3] = {0.}; -#endif -#if defined(WARPX_DIM_3D) - // Keep these double to avoid bug in single precision - double sy_new[depos_order + 3] = {0.}; - double sy_old[depos_order + 3] = {0.}; -#endif - // Keep these double to avoid bug in single precision - double sz_new[depos_order + 3] = {0.}; - double sz_old[depos_order + 3] = {0.}; - - // --- Compute shape factors - // Compute shape factors for position as they are now and at old positions - // [ijk]_new: leftmost grid point that the particle touches - const Compute_shape_factor< depos_order > compute_shape_factor; - const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; - -#if !defined(WARPX_DIM_1D_Z) const int i_new = compute_shape_factor(sx_new+1, x_new); const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new); #endif #if defined(WARPX_DIM_3D) + double sy_new[depos_order + 3] = {0.}; + double sy_old[depos_order + 3] = {0.}; const int j_new = compute_shape_factor(sy_new+1, y_new); const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new); #endif + double sz_new[depos_order + 3] = {0.}; + double sz_old[depos_order + 3] = {0.}; const int k_new = compute_shape_factor(sz_new+1, z_new); const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new); @@ -917,7 +793,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, for (int j=djl; j<=depos_order+2-dju; j++) { amrex::Real sdxi = 0._rt; for (int i=dil; i<=depos_order+1-diu; i++) { - sdxi += wqx*(sx_old[i] - sx_new[i])*( + sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*( one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k]) +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k])); amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdxi); @@ -928,7 +804,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, for (int i=dil; i<=depos_order+2-diu; i++) { amrex::Real sdyj = 0._rt; for (int j=djl; j<=depos_order+1-dju; j++) { - sdyj += wqy*(sy_old[j] - sy_new[j])*( + sdyj += wq*invdtd.y*(sy_old[j] - sy_new[j])*( one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k]) +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k])); amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdyj); @@ -939,7 +815,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, for (int i=dil; i<=depos_order+2-diu; i++) { amrex::Real sdzk = 0._rt; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k])*( + sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*( one_third*(sx_new[i]*sy_new[j] + sx_old[i]*sy_old[j]) +one_sixth*(sx_new[i]*sy_old[j] + sx_old[i]*sy_new[j])); amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdzk); @@ -952,7 +828,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, for (int k=dkl; k<=depos_order+2-dku; k++) { amrex::Real sdxi = 0._rt; for (int i=dil; i<=depos_order+1-diu; i++) { - sdxi += wqx*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]); + sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]); amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi); #if defined(WARPX_DIM_RZ) Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} @@ -973,6 +849,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k])); amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj); #if defined(WARPX_DIM_RZ) + Complex const I = Complex{0._rt, 1._rt}; Complex xy_new = xy_new0; Complex xy_mid = xy_mid0; Complex xy_old = xy_old0; @@ -980,7 +857,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode + const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xyzmin.x*dinv.x)*wq*invdtd.x/(amrex::Real)imode *(Complex(sx_new[i]*sz_new[k], 0._rt)*(xy_new - xy_mid) + Complex(sx_old[i]*sz_old[k], 0._rt)*(xy_mid - xy_old)); amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); @@ -995,7 +872,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, for (int i=dil; i<=depos_order+2-diu; i++) { Real sdzk = 0._rt; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]); + sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]); amrex::Gpu::Atomic::AddNoRet( &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} @@ -1021,7 +898,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, } amrex::Real sdzk = 0._rt; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k]); + sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k]); amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk); } #endif @@ -1047,16 +924,16 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, * \param Jx_arr,Jy_arr,Jz_arr Array4 of current density, either full array or tile. * \param np_to_deposit Number of particles for which current is deposited. * \param dt Time step for particle level - * \param dx 3D cell size + * \param dinv 3D cell size inverse * \param xyzmin Physical lower bounds of domain. * \param lo Index lower bounds of domain. * \param q species charge. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry. */ template -void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * const xp_n, - const amrex::ParticleReal * const yp_n, - const amrex::ParticleReal * const zp_n, +void doChargeConservingDepositionShapeNImplicit ([[maybe_unused]]const amrex::ParticleReal * const xp_n, + [[maybe_unused]]const amrex::ParticleReal * const yp_n, + [[maybe_unused]]const amrex::ParticleReal * const zp_n, const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, [[maybe_unused]]const amrex::ParticleReal * const uxp_n, @@ -1071,49 +948,25 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con const amrex::Array4& Jz_arr, const long np_to_deposit, const amrex::Real dt, - const std::array& dx, - const std::array xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3 lo, const amrex::Real q, - const int n_rz_azimuthal_modes) + [[maybe_unused]] const int n_rz_azimuthal_modes) { using namespace amrex; -#if !defined(WARPX_DIM_RZ) - ignore_unused(n_rz_azimuthal_modes); -#endif // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) bool const do_ionization = ion_lev; -#if !defined(WARPX_DIM_1D_Z) - Real const dxi = 1.0_rt / dx[0]; -#endif -#if !defined(WARPX_DIM_1D_Z) - Real const xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - Real const dyi = 1.0_rt / dx[1]; - Real const ymin = xyzmin[1]; -#endif - Real const dzi = 1.0_rt / dx[2]; - Real const zmin = xyzmin[2]; -#if defined(WARPX_DIM_3D) - 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) - 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]); -#elif defined(WARPX_DIM_1D_Z) - Real const invdtdz = 1.0_rt / (dt*dx[0]); - Real const invvol = 1.0_rt / (dx[2]); +#if !defined(WARPX_DIM_3D) + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; #endif -#if defined(WARPX_DIM_RZ) - Complex const I = Complex{0._rt, 1._rt}; -#endif + amrex::XDim3 const invdtd = amrex::XDim3{(1.0_rt/dt)*dinv.y*dinv.z, + (1.0_rt/dt)*dinv.x*dinv.z, + (1.0_rt/dt)*dinv.x*dinv.y}; #if !defined(WARPX_DIM_1D_Z) Real constexpr one_third = 1.0_rt / 3.0_rt; @@ -1137,7 +990,6 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con const amrex::ParticleReal gaminv = 2.0_prt/(gamma_n + gamma_np1); #endif - // wqx, wqy wqz are particle current in each direction Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; @@ -1148,24 +1000,12 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con #if !defined(WARPX_DIM_1D_Z) ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n[ip]; -#else - ignore_unused(xp_n); #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n[ip]; -#else - ignore_unused(yp_n); #endif ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n[ip]; -#if !defined(WARPX_DIM_1D_Z) - amrex::Real const wqx = wq*invdtdx; -#endif -#if defined(WARPX_DIM_3D) - amrex::Real const wqy = wq*invdtdy; -#endif - amrex::Real const wqz = wq*invdtdz; - // computes current and old position in grid units #if defined(WARPX_DIM_RZ) amrex::Real const xp_new = xp_np1; @@ -1177,51 +1017,33 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con amrex::Real const rp_new = std::sqrt(xp_new*xp_new + yp_new*yp_new); amrex::Real const rp_old = std::sqrt(xp_old*xp_old + yp_old*yp_old); amrex::Real const rp_mid = (rp_new + rp_old)/2._rt; - amrex::Real costheta_new, sintheta_new; - if (rp_new > 0._rt) { - costheta_new = xp_new/rp_new; - sintheta_new = yp_new/rp_new; - } else { - costheta_new = 1._rt; - sintheta_new = 0._rt; - } - amrex::Real costheta_mid, sintheta_mid; - if (rp_mid > 0._rt) { - costheta_mid = xp_mid/rp_mid; - sintheta_mid = yp_mid/rp_mid; - } else { - costheta_mid = 1._rt; - sintheta_mid = 0._rt; - } - amrex::Real costheta_old, sintheta_old; - if (rp_old > 0._rt) { - costheta_old = xp_old/rp_old; - sintheta_old = yp_old/rp_old; - } else { - costheta_old = 1._rt; - sintheta_old = 0._rt; - } + const amrex::Real costheta_mid = (rp_mid > 0._rt ? xp_mid/rp_mid : 1._rt); + const amrex::Real sintheta_mid = (rp_mid > 0._rt ? yp_mid/rp_mid : 0._rt); + const amrex::Real costheta_new = (rp_new > 0._rt ? xp_new/rp_new : 1._rt); + const amrex::Real sintheta_new = (rp_new > 0._rt ? yp_new/rp_new : 0._rt); + const amrex::Real costheta_old = (rp_old > 0._rt ? xp_old/rp_old : 1._rt); + const amrex::Real sintheta_old = (rp_old > 0._rt ? yp_old/rp_old : 0._rt); 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}; // Keep these double to avoid bug in single precision - double const x_new = (rp_new - xmin)*dxi; - double const x_old = (rp_old - xmin)*dxi; + double const x_new = (rp_new - xyzmin.x)*dinv.x; + double const x_old = (rp_old - xyzmin.x)*dinv.x; #else #if !defined(WARPX_DIM_1D_Z) // Keep these double to avoid bug in single precision - double const x_new = (xp_np1 - xmin)*dxi; - double const x_old = (xp_n[ip] - xmin)*dxi; + double const x_new = (xp_np1 - xyzmin.x)*dinv.x; + double const x_old = (xp_n[ip] - xyzmin.x)*dinv.x; #endif #endif #if defined(WARPX_DIM_3D) // Keep these double to avoid bug in single precision - double const y_new = (yp_np1 - ymin)*dyi; - double const y_old = (yp_n[ip] - ymin)*dyi; + double const y_new = (yp_np1 - xyzmin.y)*dinv.y; + double const y_old = (yp_n[ip] - xyzmin.y)*dinv.y; #endif // Keep these double to avoid bug in single precision - double const z_new = (zp_np1 - zmin)*dzi; - double const z_old = (zp_n[ip] - zmin)*dzi; + double const z_new = (zp_np1 - xyzmin.z)*dinv.z; + double const z_old = (zp_n[ip] - xyzmin.z)*dinv.z; #if defined(WARPX_DIM_RZ) amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv; @@ -1232,6 +1054,12 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con amrex::Real const vy = uyp_nph[ip]*gaminv; #endif + // --- Compute shape factors + // Compute shape factors for position as they are now and at old positions + // [ijk]_new: leftmost grid point that the particle touches + const Compute_shape_factor< depos_order > compute_shape_factor; + const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; + // Shape factor arrays // Note that there are extra values above and below // to possibly hold the factor for the old particle @@ -1240,30 +1068,17 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con #if !defined(WARPX_DIM_1D_Z) double sx_new[depos_order + 3] = {0.}; double sx_old[depos_order + 3] = {0.}; -#endif -#if defined(WARPX_DIM_3D) - // Keep these double to avoid bug in single precision - double sy_new[depos_order + 3] = {0.}; - double sy_old[depos_order + 3] = {0.}; -#endif - // Keep these double to avoid bug in single precision - double sz_new[depos_order + 3] = {0.}; - double sz_old[depos_order + 3] = {0.}; - - // --- Compute shape factors - // Compute shape factors for position as they are now and at old positions - // [ijk]_new: leftmost grid point that the particle touches - const Compute_shape_factor< depos_order > compute_shape_factor; - const Compute_shifted_shape_factor< depos_order > compute_shifted_shape_factor; - -#if !defined(WARPX_DIM_1D_Z) const int i_new = compute_shape_factor(sx_new+1, x_new); const int i_old = compute_shifted_shape_factor(sx_old, x_old, i_new); #endif #if defined(WARPX_DIM_3D) + double sy_new[depos_order + 3] = {0.}; + double sy_old[depos_order + 3] = {0.}; const int j_new = compute_shape_factor(sy_new+1, y_new); const int j_old = compute_shifted_shape_factor(sy_old, y_old, j_new); #endif + double sz_new[depos_order + 3] = {0.}; + double sz_old[depos_order + 3] = {0.}; const int k_new = compute_shape_factor(sz_new+1, z_new); const int k_old = compute_shifted_shape_factor(sz_old, z_old, k_new); @@ -1288,7 +1103,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con for (int j=djl; j<=depos_order+2-dju; j++) { amrex::Real sdxi = 0._rt; for (int i=dil; i<=depos_order+1-diu; i++) { - sdxi += wqx*(sx_old[i] - sx_new[i])*( + sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*( one_third*(sy_new[j]*sz_new[k] + sy_old[j]*sz_old[k]) +one_sixth*(sy_new[j]*sz_old[k] + sy_old[j]*sz_new[k])); amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdxi); @@ -1299,7 +1114,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con for (int i=dil; i<=depos_order+2-diu; i++) { amrex::Real sdyj = 0._rt; for (int j=djl; j<=depos_order+1-dju; j++) { - sdyj += wqy*(sy_old[j] - sy_new[j])*( + sdyj += wq*invdtd.y*(sy_old[j] - sy_new[j])*( one_third*(sx_new[i]*sz_new[k] + sx_old[i]*sz_old[k]) +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k])); amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdyj); @@ -1310,7 +1125,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con for (int i=dil; i<=depos_order+2-diu; i++) { amrex::Real sdzk = 0._rt; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k])*( + sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*( one_third*(sx_new[i]*sy_new[j] + sx_old[i]*sy_old[j]) +one_sixth*(sx_new[i]*sy_old[j] + sx_old[i]*sy_new[j])); amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+i_new-1+i, lo.y+j_new-1+j, lo.z+k_new-1+k), sdzk); @@ -1323,7 +1138,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con for (int k=dkl; k<=depos_order+2-dku; k++) { amrex::Real sdxi = 0._rt; for (int i=dil; i<=depos_order+1-diu; i++) { - sdxi += wqx*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]); + sdxi += wq*invdtd.x*(sx_old[i] - sx_new[i])*0.5_rt*(sz_new[k] + sz_old[k]); amrex::Gpu::Atomic::AddNoRet( &Jx_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdxi); #if defined(WARPX_DIM_RZ) Complex xy_mid = xy_mid0; // Throughout the following loop, xy_mid takes the value e^{i m theta} @@ -1344,6 +1159,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con +one_sixth*(sx_new[i]*sz_old[k] + sx_old[i]*sz_new[k])); amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 0), sdyj); #if defined(WARPX_DIM_RZ) + Complex const I = Complex{0._rt, 1._rt}; Complex xy_new = xy_new0; Complex xy_mid = xy_mid0; Complex xy_old = xy_old0; @@ -1351,7 +1167,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con for (int imode=1 ; imode < n_rz_azimuthal_modes ; imode++) { // The factor 2 comes from the normalization of the modes // The minus sign comes from the different convention with respect to Davidson et al. - const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xmin*dxi)*wq*invdtdx/(amrex::Real)imode + const Complex djt_cmplx = -2._rt * I*(i_new-1 + i + xyzmin.x*dinv.x)*wq*invdtd.x/(amrex::Real)imode *(Complex(sx_new[i]*sz_new[k], 0._rt)*(xy_new - xy_mid) + Complex(sx_old[i]*sz_old[k], 0._rt)*(xy_mid - xy_old)); amrex::Gpu::Atomic::AddNoRet( &Jy_arr(lo.x+i_new-1+i, lo.y+k_new-1+k, 0, 2*imode-1), djt_cmplx.real()); @@ -1366,7 +1182,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con for (int i=dil; i<=depos_order+2-diu; i++) { Real sdzk = 0._rt; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]); + sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k])*0.5_rt*(sx_new[i] + sx_old[i]); amrex::Gpu::Atomic::AddNoRet( &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} @@ -1392,7 +1208,7 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con } amrex::Real sdzk = 0._rt; for (int k=dkl; k<=depos_order+1-dku; k++) { - sdzk += wqz*(sz_old[k] - sz_new[k]); + sdzk += wq*invdtd.z*(sz_old[k] - sz_new[k]); amrex::Gpu::Atomic::AddNoRet( &Jz_arr(lo.x+k_new-1+k, 0, 0, 0), sdzk); } #endif @@ -1420,16 +1236,16 @@ void doChargeConservingDepositionShapeNImplicit (const amrex::ParticleReal * con * \param Jx_arr,Jy_arr,Jz_arr Array4 of current density, either full array or tile. * \param np_to_deposit Number of particles for which current is deposited. * \param dt Time step for particle level - * \param dx 3D cell size + * \param dinv 3D cell size inverse * \param xyzmin Physical lower bounds of domain. * \param lo Index lower bounds of domain. * \param q species charge. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry. */ template -void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_n, - const amrex::ParticleReal * const yp_n, - const amrex::ParticleReal * const zp_n, +void doVillasenorDepositionShapeNImplicit ([[maybe_unused]]const amrex::ParticleReal * const xp_n, + [[maybe_unused]]const amrex::ParticleReal * const yp_n, + [[maybe_unused]]const amrex::ParticleReal * const zp_n, const GetParticlePosition& GetPosition, const amrex::ParticleReal * const wp, [[maybe_unused]]const amrex::ParticleReal * const uxp_n, @@ -1444,42 +1260,21 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ const amrex::Array4& Jz_arr, const long np_to_deposit, const amrex::Real dt, - const std::array& dx, - const std::array xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3 lo, const amrex::Real q, - const int n_rz_azimuthal_modes) + [[maybe_unused]] const int n_rz_azimuthal_modes) { using namespace amrex; -#if !defined(WARPX_DIM_RZ) - ignore_unused(n_rz_azimuthal_modes); -#endif // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer // (do_ionization=1) bool const do_ionization = ion_lev; -#if !defined(WARPX_DIM_1D_Z) - Real const dxi = 1.0_rt / dx[0]; -#endif -#if !defined(WARPX_DIM_1D_Z) - Real const xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - Real const dyi = 1.0_rt / dx[1]; - Real const ymin = xyzmin[1]; -#endif - Real const dzi = 1.0_rt / dx[2]; - Real const zmin = xyzmin[2]; -#if defined(WARPX_DIM_3D) - Real const invvol = 1.0_rt / (dx[0]*dx[1]*dx[2]); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - Real const invvol = 1.0_rt / (dx[0]*dx[2]); -#elif defined(WARPX_DIM_1D_Z) - Real const invvol = 1.0_rt / (dx[2]); -#endif + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; -#if !defined(WARPX_DIM_1D_Z) +#if (AMREX_SPACEDIM > 1) Real constexpr one_third = 1.0_rt / 3.0_rt; Real constexpr one_sixth = 1.0_rt / 6.0_rt; #endif @@ -1502,7 +1297,6 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ const amrex::ParticleReal gaminv = 2.0_prt/(gamma_n + gamma_np1); #endif - // wqx, wqy wqz are particle current in each direction Real wq = q*wp[ip]; if (do_ionization){ wq *= ion_lev[ip]; @@ -1513,13 +1307,9 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ #if !defined(WARPX_DIM_1D_Z) ParticleReal const xp_np1 = 2._prt*xp_nph - xp_n[ip]; -#else - ignore_unused(xp_n); #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) ParticleReal const yp_np1 = 2._prt*yp_nph - yp_n[ip]; -#else - ignore_unused(yp_n); #endif ParticleReal const zp_np1 = 2._prt*zp_nph - zp_n[ip]; @@ -1545,14 +1335,14 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid}; // Keep these double to avoid bug in single precision - double const x_new = (rp_new - xmin)*dxi; - double const x_old = (rp_old - xmin)*dxi; + double const x_new = (rp_new - xyzmin.x)*dinv.x; + double const x_old = (rp_old - xyzmin.x)*dinv.x; amrex::Real const vx = (rp_new - rp_old)/dt; amrex::Real const vy = (-uxp_nph[ip]*sintheta_mid + uyp_nph[ip]*costheta_mid)*gaminv; #elif defined(WARPX_DIM_XZ) // Keep these double to avoid bug in single precision - double const x_new = (xp_np1 - xmin)*dxi; - double const x_old = (xp_n[ip] - xmin)*dxi; + double const x_new = (xp_np1 - xyzmin.x)*dinv.x; + double const x_old = (xp_n[ip] - xyzmin.x)*dinv.x; amrex::Real const vx = (xp_np1 - xp_n[ip])/dt; amrex::Real const vy = uyp_nph[ip]*gaminv; #elif defined(WARPX_DIM_1D_Z) @@ -1560,17 +1350,17 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ amrex::Real const vy = uyp_nph[ip]*gaminv; #elif defined(WARPX_DIM_3D) // Keep these double to avoid bug in single precision - double const x_new = (xp_np1 - xmin)*dxi; - double const x_old = (xp_n[ip] - xmin)*dxi; - double const y_new = (yp_np1 - ymin)*dyi; - double const y_old = (yp_n[ip] - ymin)*dyi; + double const x_new = (xp_np1 - xyzmin.x)*dinv.x; + double const x_old = (xp_n[ip] - xyzmin.x)*dinv.x; + double const y_new = (yp_np1 - xyzmin.y)*dinv.y; + double const y_old = (yp_n[ip] - xyzmin.y)*dinv.y; amrex::Real const vx = (xp_np1 - xp_n[ip])/dt; amrex::Real const vy = (yp_np1 - yp_n[ip])/dt; #endif // Keep these double to avoid bug in single precision - double const z_new = (zp_np1 - zmin)*dzi; - double const z_old = (zp_n[ip] - zmin)*dzi; + double const z_new = (zp_np1 - xyzmin.z)*dinv.z; + double const z_old = (zp_n[ip] - xyzmin.z)*dinv.z; amrex::Real const vz = (zp_np1 - zp_n[ip])/dt; // Define velocity kernals to deposit @@ -2057,7 +1847,7 @@ void doVillasenorDepositionShapeNImplicit (const amrex::ParticleReal * const xp_ * current positions of the particles. When different than 0, * the particle position will be temporarily modified to match * the time of the deposition. - * \param[in] dx 3D cell size + * \param[in] dinv 3D cell size inverse * \param[in] xyzmin 3D lower bounds of physical domain * \param[in] lo Dimension-agnostic lower bounds of index domain * \param[in] q Species charge @@ -2076,57 +1866,37 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, long np_to_deposit, amrex::Real dt, amrex::Real relative_time, - const std::array& dx, - const std::array& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, amrex::Dim3 lo, amrex::Real q, - int n_rz_azimuthal_modes) + [[maybe_unused]]int n_rz_azimuthal_modes) { using namespace amrex::literals; #if defined(WARPX_DIM_RZ) amrex::ignore_unused(GetPosition, wp, uxp, uyp, uzp, ion_lev, Dx_fab, Dy_fab, Dz_fab, - np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes); + np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q); WARPX_ABORT_WITH_MESSAGE("Vay deposition not implemented in RZ geometry"); #endif #if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(GetPosition, wp, uxp, uyp, uzp, ion_lev, Dx_fab, Dy_fab, Dz_fab, - np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, n_rz_azimuthal_modes); - WARPX_ABORT_WITH_MESSAGE("Vay deposition not implemented in cartesian 1D geometry"); + np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q); + WARPX_ABORT_WITH_MESSAGE("Vay deposition not implemented in 1D geometry"); #endif #if !(defined WARPX_DIM_RZ || defined WARPX_DIM_1D_Z) - amrex::ignore_unused(n_rz_azimuthal_modes); // If ion_lev is a null pointer, then do_ionization=0, else do_ionization=1 const bool do_ionization = ion_lev; - // Inverse cell volume in each direction - const amrex::Real dxi = 1._rt / dx[0]; - const amrex::Real dzi = 1._rt / dx[2]; -#if defined(WARPX_DIM_3D) - const amrex::Real dyi = 1._rt / dx[1]; -#endif - // Inverse of time step const amrex::Real invdt = 1._rt / dt; - // Total inverse cell volume -#if defined(WARPX_DIM_XZ) - const amrex::Real invvol = dxi * dzi; -#elif defined(WARPX_DIM_3D) - const amrex::Real invvol = dxi * dyi * dzi; -#endif - - // Lower bound of physical domain in each direction - const amrex::Real xmin = xyzmin[0]; - const amrex::Real zmin = xyzmin[2]; -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; // Allocate temporary arrays #if defined(WARPX_DIM_3D) @@ -2172,23 +1942,18 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, yp += relative_time * vy; zp += relative_time * vz; - // Particle current densities -#if defined(WARPX_DIM_XZ) - const amrex::Real wqy = wq * vy * invvol; -#endif - // Current and old particle positions in grid units // Keep these double to avoid bug in single precision. - double const x_new = (xp - xmin + 0.5_rt*dt*vx) * dxi; - double const x_old = (xp - xmin - 0.5_rt*dt*vx) * dxi; + double const x_new = (xp - xyzmin.x + 0.5_rt*dt*vx) * dinv.x; + double const x_old = (xp - xyzmin.x - 0.5_rt*dt*vx) * dinv.x; #if defined(WARPX_DIM_3D) // Keep these double to avoid bug in single precision. - double const y_new = (yp - ymin + 0.5_rt*dt*vy) * dyi; - double const y_old = (yp - ymin - 0.5_rt*dt*vy) * dyi; + double const y_new = (yp - xyzmin.y + 0.5_rt*dt*vy) * dinv.y; + double const y_old = (yp - xyzmin.y - 0.5_rt*dt*vy) * dinv.y; #endif // Keep these double to avoid bug in single precision. - double const z_new = (zp - zmin + 0.5_rt*dt*vz) * dzi; - double const z_old = (zp - zmin - 0.5_rt*dt*vz) * dzi; + double const z_new = (zp - xyzmin.z + 0.5_rt*dt*vz) * dinv.z; + double const z_old = (zp - xyzmin.z - 0.5_rt*dt*vz) * dinv.z; // Shape factor arrays for current and old positions (nodal) // Keep these double to avoid bug in single precision. @@ -2235,6 +2000,7 @@ void doVayDepositionShapeN (const GetParticlePosition& GetPosition, // Deposit current into Dx_arr, Dy_arr and Dz_arr #if defined(WARPX_DIM_XZ) + const amrex::Real wqy = wq * vy * invvol; for (int k=0; k<=depos_order; k++) { for (int i=0; i<=depos_order; i++) { diff --git a/Source/Particles/Deposition/SharedDepositionUtils.H b/Source/Particles/Deposition/SharedDepositionUtils.H index 1e2294be3a2..e6deb51b5df 100644 --- a/Source/Particles/Deposition/SharedDepositionUtils.H +++ b/Source/Particles/Deposition/SharedDepositionUtils.H @@ -74,8 +74,8 @@ void depositComponent (const GetParticlePosition& GetPosition, amrex::Array4 const& j_buff, amrex::IntVect const j_type, const amrex::Real relative_time, - const std::array& dx, - const std::array& xyzmin, + const amrex::XDim3 dinv, + const amrex::XDim3 xyzmin, const amrex::Dim3 lo, const amrex::Real q, const int n_rz_azimuthal_modes, @@ -91,26 +91,8 @@ void depositComponent (const GetParticlePosition& GetPosition, // 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 dzi = 1.0_rt/dx[2]; -#if defined(WARPX_DIM_1D_Z) - const amrex::Real invvol = dzi; -#endif -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real invvol = dxi*dzi; -#elif defined(WARPX_DIM_3D) - const amrex::Real dxi = 1.0_rt/dx[0]; - const amrex::Real dyi = 1.0_rt/dx[1]; - const amrex::Real invvol = dxi*dyi*dzi; -#endif -#if (AMREX_SPACEDIM >= 2) - const amrex::Real xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif - const amrex::Real zmin = xyzmin[2]; + const amrex::Real invvol = dinv.x*dinv.y*dinv.z; const amrex::Real clightsq = 1.0_rt/PhysConst::c/PhysConst::c; @@ -171,9 +153,9 @@ void depositComponent (const GetParticlePosition& GetPosition, // Get particle position after 1/2 push back in position #if defined(WARPX_DIM_RZ) // Keep these double to avoid bug in single precision - const double xmid = (rpmid - xmin)*dxi; + const double xmid = (rpmid - xyzmin.x)*dinv.x; #else - const double xmid = ((xp - xmin) + relative_time*vx)*dxi; + const double xmid = ((xp - xyzmin.x) + relative_time*vx)*dinv.x; #endif // j_j[xyz] leftmost grid point in x that the particle touches for the centering of each current // sx_j[xyz] shape factor along x for the centering of each current @@ -203,7 +185,7 @@ void depositComponent (const GetParticlePosition& GetPosition, #if defined(WARPX_DIM_3D) // y direction // Keep these double to avoid bug in single precision - const double ymid = ((yp - ymin) + relative_time*vy)*dyi; + const double ymid = ((yp - xyzmin.y) + relative_time*vy)*dinv.y; double sy_node[depos_order + 1] = {0.}; double sy_cell[depos_order + 1] = {0.}; int k_node = 0; @@ -224,7 +206,7 @@ void depositComponent (const GetParticlePosition& GetPosition, // z direction // Keep these double to avoid bug in single precision - const double zmid = ((zp - zmin) + relative_time*vz)*dzi; + const double zmid = ((zp - xyzmin.z) + relative_time*vz)*dinv.z; double sz_node[depos_order + 1] = {0.}; double sz_cell[depos_order + 1] = {0.}; int l_node = 0; diff --git a/Source/Particles/ElementaryProcess/Ionization.H b/Source/Particles/ElementaryProcess/Ionization.H index cee7ec07eb5..6f98c18959a 100644 --- a/Source/Particles/ElementaryProcess/Ionization.H +++ b/Source/Particles/ElementaryProcess/Ionization.H @@ -62,8 +62,8 @@ struct IonizationFilterFunc amrex::IndexType m_by_type; amrex::IndexType m_bz_type; - amrex::GpuArray m_dx_arr; - amrex::GpuArray m_xyzmin_arr; + amrex::XDim3 m_dinv; + amrex::XDim3 m_xyzmin; bool m_galerkin_interpolation; int m_nox; @@ -117,7 +117,7 @@ struct IonizationFilterFunc doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz, m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr, m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type, - m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes, + m_dinv, m_xyzmin, m_lo, m_n_rz_azimuthal_modes, m_nox, m_galerkin_interpolation); m_get_externalEB(i, ex, ey, ez, bx, by, bz); diff --git a/Source/Particles/ElementaryProcess/Ionization.cpp b/Source/Particles/ElementaryProcess/Ionization.cpp index b7b91e4d4e3..0568e302eec 100644 --- a/Source/Particles/ElementaryProcess/Ionization.cpp +++ b/Source/Particles/ElementaryProcess/Ionization.cpp @@ -76,11 +76,10 @@ IonizationFilterFunc::IonizationFilterFunc (const WarpXParIter& a_pti, int lev, box.grow(ngEB); const std::array& dx = WarpX::CellSize(std::max(lev, 0)); - m_dx_arr = {dx[0], dx[1], dx[2]}; + m_dinv = amrex::XDim3{1._rt/dx[0], 1._rt/dx[1], 1._rt/dx[2]}; // Lower corner of tile box physical domain (take into account Galilean shift) - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + m_xyzmin = WarpX::LowerCorner(box, lev, 0._rt); m_lo = amrex::lbound(box); } diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.H b/Source/Particles/ElementaryProcess/QEDPairGeneration.H index d550c5f03ee..180fdf0fb35 100644 --- a/Source/Particles/ElementaryProcess/QEDPairGeneration.H +++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.H @@ -144,7 +144,7 @@ public: doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz, m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr, m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type, - m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes, + m_dinv, m_xyzmin, m_lo, m_n_rz_azimuthal_modes, m_nox, m_galerkin_interpolation); //Despite the names of the variables, positrons and electrons @@ -198,8 +198,8 @@ private: amrex::IndexType m_by_type; amrex::IndexType m_bz_type; - amrex::GpuArray m_dx_arr; - amrex::GpuArray m_xyzmin_arr; + amrex::XDim3 m_dinv; + amrex::XDim3 m_xyzmin; bool m_galerkin_interpolation; int m_nox; diff --git a/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp b/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp index 2b380d454f4..82546e43ef4 100644 --- a/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp +++ b/Source/Particles/ElementaryProcess/QEDPairGeneration.cpp @@ -64,11 +64,10 @@ PairGenerationTransformFunc (BreitWheelerGeneratePairs const generate_functor, box.grow(ngEB); const std::array& dx = WarpX::CellSize(std::max(lev, 0)); - m_dx_arr = {dx[0], dx[1], dx[2]}; + m_dinv = amrex::XDim3{1._rt/dx[0], 1._rt/dx[1], 1._rt/dx[2]}; // Lower corner of tile box physical domain (take into account Galilean shift) - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + m_xyzmin = WarpX::LowerCorner(box, lev, 0._rt); m_lo = amrex::lbound(box); } diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H index dd5f57a8f24..514526374bd 100644 --- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.H +++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.H @@ -156,7 +156,7 @@ public: doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz, m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr, m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type, - m_dx_arr, m_xyzmin_arr, m_lo, m_n_rz_azimuthal_modes, + m_dinv, m_xyzmin, m_lo, m_n_rz_azimuthal_modes, m_nox, m_galerkin_interpolation); auto& ux = src.m_rdata[PIdx::ux][i_src]; @@ -209,8 +209,8 @@ private: amrex::IndexType m_by_type; amrex::IndexType m_bz_type; - amrex::GpuArray m_dx_arr; - amrex::GpuArray m_xyzmin_arr; + amrex::XDim3 m_dinv; + amrex::XDim3 m_xyzmin; bool m_galerkin_interpolation; int m_nox; diff --git a/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp b/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp index 077a4659ce5..688c3c0184d 100644 --- a/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp +++ b/Source/Particles/ElementaryProcess/QEDPhotonEmission.cpp @@ -67,13 +67,10 @@ PhotonEmissionTransformFunc (QuantumSynchrotronGetOpticalDepth opt_depth_functor box.grow(ngEB); const std::array& dx = WarpX::CellSize(std::max(lev, 0)); - m_dx_arr = {dx[0], dx[1], dx[2]}; + m_dinv = amrex::XDim3{1._rt/dx[0], 1._rt/dx[1], 1._rt/dx[2]}; // Lower corner of tile box physical domain (take into account Galilean shift) - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); - m_xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - - + m_xyzmin = WarpX::LowerCorner(box, lev, 0._rt); m_lo = amrex::lbound(box); } diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index b5bd4376ba1..4b4590b8642 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -28,16 +28,16 @@ * \param bx_arr,by_arr,bz_arr Array4 of the magnetic field, either full array or tile. * \param ex_type,ey_type,ez_type IndexType of the electric field * \param bx_type,by_type,bz_type IndexType of the magnetic field - * \param dx 3D cell spacing - * \param xyzmin Physical lower bounds of domain in x, y, z. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry */ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -void doGatherShapeN (const amrex::ParticleReal xp, - const amrex::ParticleReal yp, - const amrex::ParticleReal zp, +void doGatherShapeN ([[maybe_unused]] const amrex::ParticleReal xp, + [[maybe_unused]] const amrex::ParticleReal yp, + [[maybe_unused]] const amrex::ParticleReal zp, amrex::ParticleReal& Exp, amrex::ParticleReal& Eyp, amrex::ParticleReal& Ezp, @@ -56,41 +56,13 @@ void doGatherShapeN (const amrex::ParticleReal xp, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, - const amrex::GpuArray& dx, - const amrex::GpuArray& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3& lo, - const int n_rz_azimuthal_modes) + [[maybe_unused]] const int n_rz_azimuthal_modes) { using namespace amrex; -#if defined(WARPX_DIM_XZ) - amrex::ignore_unused(yp); -#endif - -#if defined(WARPX_DIM_1D_Z) - amrex::ignore_unused(xp,yp); -#endif - -#ifndef WARPX_DIM_RZ - amrex::ignore_unused(n_rz_azimuthal_modes); -#endif - -#if (AMREX_SPACEDIM >= 2) - const amrex::Real dxi = 1.0_rt/dx[0]; -#endif - const amrex::Real dzi = 1.0_rt/dx[2]; -#if defined(WARPX_DIM_3D) - const amrex::Real dyi = 1.0_rt/dx[1]; -#endif - -#if (AMREX_SPACEDIM >= 2) - const amrex::Real xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - const amrex::Real ymin = xyzmin[1]; -#endif - const amrex::Real zmin = xyzmin[2]; - constexpr int zdir = WARPX_ZINDEX; constexpr int NODE = amrex::IndexType::NODE; constexpr int CELL = amrex::IndexType::CELL; @@ -105,9 +77,9 @@ void doGatherShapeN (const amrex::ParticleReal xp, // Get particle position #ifdef WARPX_DIM_RZ const amrex::Real rp = std::sqrt(xp*xp + yp*yp); - const amrex::Real x = (rp - xmin)*dxi; + const amrex::Real x = (rp - xyzmin.x)*dinv.x; #else - const amrex::Real x = (xp-xmin)*dxi; + const amrex::Real x = (xp-xyzmin.x)*dinv.x; #endif // j_[eb][xyz] leftmost grid point in x that the particle touches for the centering of each current @@ -151,7 +123,7 @@ void doGatherShapeN (const amrex::ParticleReal xp, #if defined(WARPX_DIM_3D) // y direction - const amrex::Real y = (yp-ymin)*dyi; + const amrex::Real y = (yp-xyzmin.y)*dinv.y; amrex::Real sy_node[depos_order + 1]; amrex::Real sy_cell[depos_order + 1]; amrex::Real sy_node_v[depos_order + 1 - galerkin_interpolation]; @@ -187,7 +159,7 @@ void doGatherShapeN (const amrex::ParticleReal xp, #endif // z direction - const amrex::Real z = (zp-zmin)*dzi; + const amrex::Real z = (zp-xyzmin.z)*dinv.z; amrex::Real sz_node[depos_order + 1]; amrex::Real sz_cell[depos_order + 1]; amrex::Real sz_node_v[depos_order + 1 - galerkin_interpolation]; @@ -464,10 +436,10 @@ void doGatherShapeN (const amrex::ParticleReal xp, * \param Bx_arr,By_arr,Bz_arr Array4 of the magnetic field, either full array or tile. * \param Ex_type,Ey_type,Ez_type IndexType of the electric field * \param Bx_type,By_type,Bz_type IndexType of the magnetic field - * \param dx 3D cell spacing - * \param xyzmin Physical lower bounds of domain in x, y, z. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. - * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry + * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry */ template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE @@ -496,8 +468,8 @@ void doGatherShapeNEsirkepovStencilImplicit ( [[maybe_unused]] const amrex::IndexType Bx_type, [[maybe_unused]] const amrex::IndexType By_type, [[maybe_unused]] const amrex::IndexType Bz_type, - const amrex::GpuArray& dx, - const amrex::GpuArray& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3& lo, const int n_rz_azimuthal_modes) { @@ -506,19 +478,6 @@ void doGatherShapeNEsirkepovStencilImplicit ( ignore_unused(n_rz_azimuthal_modes); #endif -#if !defined(WARPX_DIM_1D_Z) - Real const dxi = 1.0_rt / dx[0]; -#endif -#if !defined(WARPX_DIM_1D_Z) - Real const xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - Real const dyi = 1.0_rt / dx[1]; - Real const ymin = xyzmin[1]; -#endif - Real const dzi = 1.0_rt / dx[2]; - Real const zmin = xyzmin[2]; - #if !defined(WARPX_DIM_1D_Z) Real constexpr one_third = 1.0_rt / 3.0_rt; Real constexpr one_sixth = 1.0_rt / 6.0_rt; @@ -553,23 +512,23 @@ void doGatherShapeNEsirkepovStencilImplicit ( } const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid}; // Keep these double to avoid bug in single precision - double const x_new = (rp_new - xmin)*dxi; - double const x_old = (rp_old - xmin)*dxi; + double const x_new = (rp_new - xyzmin.x)*dinv.x; + double const x_old = (rp_old - xyzmin.x)*dinv.x; #else #if !defined(WARPX_DIM_1D_Z) // Keep these double to avoid bug in single precision - double const x_new = (xp_np1 - xmin)*dxi; - double const x_old = (xp_n - xmin)*dxi; + double const x_new = (xp_np1 - xyzmin.x)*dinv.x; + double const x_old = (xp_n - xyzmin.x)*dinv.x; #endif #endif #if defined(WARPX_DIM_3D) // Keep these double to avoid bug in single precision - double const y_new = (yp_np1 - ymin)*dyi; - double const y_old = (yp_n - ymin)*dyi; + double const y_new = (yp_np1 - xyzmin.y)*dinv.y; + double const y_old = (yp_n - xyzmin.y)*dinv.y; #endif // Keep these double to avoid bug in single precision - double const z_new = (zp_np1 - zmin)*dzi; - double const z_old = (zp_n - zmin)*dzi; + double const z_new = (zp_np1 - xyzmin.z)*dinv.z; + double const z_old = (zp_n - xyzmin.z)*dinv.z; // Shape factor arrays // Note that there are extra values above and below @@ -894,8 +853,8 @@ void doGatherShapeNEsirkepovStencilImplicit ( * \param Bx_arr,By_arr,Bz_arr Array4 of the magnetic field, either full array or tile. * \param Ex_type,Ey_type,Ez_type IndexType of the electric field * \param Bx_type,By_type,Bz_type IndexType of the magnetic field - * \param dx 3D cell spacing - * \param xyzmin Physical lower bounds of domain in x, y, z. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry */ @@ -926,8 +885,8 @@ void doGatherPicnicShapeN ( [[maybe_unused]] const amrex::IndexType Bx_type, [[maybe_unused]] const amrex::IndexType By_type, [[maybe_unused]] const amrex::IndexType Bz_type, - const amrex::GpuArray& dx, - const amrex::GpuArray& xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3& lo, const int n_rz_azimuthal_modes) { @@ -936,19 +895,6 @@ void doGatherPicnicShapeN ( ignore_unused(n_rz_azimuthal_modes); #endif -#if !defined(WARPX_DIM_1D_Z) - Real const dxi = 1.0_rt / dx[0]; -#endif -#if !defined(WARPX_DIM_1D_Z) - Real const xmin = xyzmin[0]; -#endif -#if defined(WARPX_DIM_3D) - Real const dyi = 1.0_rt / dx[1]; - Real const ymin = xyzmin[1]; -#endif - Real const dzi = 1.0_rt / dx[2]; - Real const zmin = xyzmin[2]; - #if !defined(WARPX_DIM_1D_Z) const ParticleReal xp_np1 = 2._prt*xp_nph - xp_n; #endif @@ -983,25 +929,25 @@ void doGatherPicnicShapeN ( } const Complex xy_mid0 = Complex{costheta_mid, sintheta_mid}; // Keep these double to avoid bug in single precision - double const x_new = (rp_new - xmin)*dxi; - double const x_old = (rp_old - xmin)*dxi; - double const x_bar = (rp_mid - xmin)*dxi; + double const x_new = (rp_new - xyzmin.x)*dinv.x; + double const x_old = (rp_old - xyzmin.x)*dinv.x; + double const x_bar = (rp_mid - xyzmin.x)*dinv.x; #elif !defined(WARPX_DIM_1D_Z) // Keep these double to avoid bug in single precision - double const x_new = (xp_np1 - xmin)*dxi; - double const x_old = (xp_n - xmin)*dxi; - double const x_bar = (xp_nph - xmin)*dxi; + double const x_new = (xp_np1 - xyzmin.x)*dinv.x; + double const x_old = (xp_n - xyzmin.x)*dinv.x; + double const x_bar = (xp_nph - xyzmin.x)*dinv.x; #endif #if defined(WARPX_DIM_3D) // Keep these double to avoid bug in single precision - double const y_new = (yp_np1 - ymin)*dyi; - double const y_old = (yp_n - ymin)*dyi; - double const y_bar = (yp_nph - ymin)*dyi; + double const y_new = (yp_np1 - xyzmin.y)*dinv.y; + double const y_old = (yp_n - xyzmin.y)*dinv.y; + double const y_bar = (yp_nph - xyzmin.y)*dinv.y; #endif // Keep these double to avoid bug in single precision - double const z_new = (zp_np1 - zmin)*dzi; - double const z_old = (zp_n - zmin)*dzi; - double const z_bar = (zp_nph - zmin)*dzi; + double const z_new = (zp_np1 - xyzmin.z)*dinv.z; + double const z_old = (zp_n - xyzmin.z)*dinv.z; + double const z_bar = (zp_nph - xyzmin.z)*dinv.z; // 1) Determine the number of segments. // 2) Loop over segments and gather electric field. @@ -1567,8 +1513,8 @@ void doGatherPicnicShapeN ( * \param exfab,eyfab,ezfab Array4 of the electric field, either full array or tile. * \param bxfab,byfab,bzfab Array4 of the magnetic field, either full array or tile. * \param np_to_gather Number of particles for which field is gathered. - * \param dx 3D cell size - * \param xyzmin Physical lower bounds of domain. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry */ @@ -1585,15 +1531,12 @@ void doGatherShapeN(const GetParticlePosition& getPosition, amrex::FArrayBox const * const byfab, amrex::FArrayBox const * const bzfab, const long np_to_gather, - const std::array& dx, - const std::array xyzmin, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3 lo, const int n_rz_azimuthal_modes) { - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); amrex::Array4 const& ez_arr = ezfab->array(); @@ -1622,7 +1565,7 @@ void doGatherShapeN(const GetParticlePosition& getPosition, xp, yp, zp, Exp[ip], Eyp[ip], Ezp[ip], Bxp[ip], Byp[ip], Bzp[ip], ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } ); } @@ -1637,8 +1580,8 @@ void doGatherShapeN(const GetParticlePosition& getPosition, * \param bx_arr,by_arr,bz_arr Array4 of the magnetic field, either full array or tile. * \param ex_type,ey_type,ez_type IndexType of the electric field * \param bx_type,by_type,bz_type IndexType of the magnetic field - * \param dx_arr 3D cell spacing - * \param xyzmin_arr Physical lower bounds of domain in x, y, z. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry * \param nox order of the particle shape function @@ -1666,8 +1609,8 @@ void doGatherShapeN (const amrex::ParticleReal xp, const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, - const amrex::GpuArray& dx_arr, - const amrex::GpuArray& xyzmin_arr, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3& lo, const int n_rz_azimuthal_modes, const int nox, @@ -1678,44 +1621,44 @@ void doGatherShapeN (const amrex::ParticleReal xp, doGatherShapeN<1,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 2) { doGatherShapeN<2,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 3) { doGatherShapeN<3,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 4) { doGatherShapeN<4,1>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } } else { if (nox == 1) { doGatherShapeN<1,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 2) { doGatherShapeN<2,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 3) { doGatherShapeN<3,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 4) { doGatherShapeN<4,0>(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } } } @@ -1732,8 +1675,8 @@ void doGatherShapeN (const amrex::ParticleReal xp, * \param bx_arr,by_arr,bz_arr Array4 of the magnetic field, either full array or tile. * \param ex_type,ey_type,ez_type IndexType of the electric field * \param bx_type,by_type,bz_type IndexType of the magnetic field - * \param dx_arr 3D cell spacing - * \param xyzmin_arr Physical lower bounds of domain in x, y, z. + * \param dinv 3D cell size inverse + * \param xyzmin The lower bounds of the domain * \param lo Index lower bounds of domain. * \param n_rz_azimuthal_modes Number of azimuthal modes when using RZ geometry * \param nox order of the particle shape function @@ -1766,8 +1709,8 @@ void doGatherShapeNImplicit ( const amrex::IndexType bx_type, const amrex::IndexType by_type, const amrex::IndexType bz_type, - const amrex::GpuArray& dx_arr, - const amrex::GpuArray& xyzmin_arr, + const amrex::XDim3 & dinv, + const amrex::XDim3 & xyzmin, const amrex::Dim3& lo, const int n_rz_azimuthal_modes, const int nox, @@ -1779,25 +1722,25 @@ void doGatherShapeNImplicit ( Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 2) { doGatherShapeNEsirkepovStencilImplicit<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 3) { doGatherShapeNEsirkepovStencilImplicit<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 4) { doGatherShapeNEsirkepovStencilImplicit<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } } else if (depos_type==3) { // CurrentDepositionAlgo::Villasenor @@ -1806,25 +1749,25 @@ void doGatherShapeNImplicit ( Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 2) { doGatherPicnicShapeN<2>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 3) { doGatherPicnicShapeN<3>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 4) { doGatherPicnicShapeN<4>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } } else if (depos_type==1) { // CurrentDepositionAlgo::Direct @@ -1832,22 +1775,22 @@ void doGatherShapeNImplicit ( doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 2) { doGatherShapeN<2,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 3) { doGatherShapeN<3,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } else if (nox == 4) { doGatherShapeN<4,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes); + dinv, xyzmin, lo, n_rz_azimuthal_modes); } } } diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index aa9f04be224..1f15d5210f5 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -93,8 +93,8 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti, int lev, int gather_lev, amrex::Real dt, ScaleFields /*scaleFields*/, DtType a_dt_type) { - // Get cell size on gather_lev - const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); + // Get inverse cell size on gather_lev + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(gather_lev,0)); // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. @@ -142,7 +142,7 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti, const amrex::ParticleReal Bz_external_particle = m_B_external_particle[2]; // Lower corner of tile box physical domain (take into account Galilean shift) - const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt); const Dim3 lo = lbound(box); @@ -150,9 +150,6 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti, const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); amrex::Array4 const& ez_arr = ezfab->array(); @@ -201,7 +198,7 @@ PhotonParticleContainer::PushPX (WarpXParIter& pti, doGatherShapeN(x, y, z, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, + dinv, xyzmin, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index e2be4f948ca..89ae435882b 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2558,7 +2558,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, if (do_not_push) { return; } - const std::array& dx = WarpX::CellSize(std::max(lev,0)); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(lev,0)); #ifdef AMREX_USE_OMP #pragma omp parallel @@ -2590,7 +2590,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const amrex::ParticleReal By_external_particle = m_B_external_particle[1]; const amrex::ParticleReal Bz_external_particle = m_B_external_particle[2]; - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt); const Dim3 lo = lbound(box); @@ -2598,9 +2598,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - amrex::Array4 const& ex_arr = exfab.array(); amrex::Array4 const& ey_arr = eyfab.array(); amrex::Array4 const& ez_arr = ezfab.array(); @@ -2657,7 +2654,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, + dinv, xyzmin, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); } @@ -2755,7 +2752,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, if (np_to_push == 0) { return; } // Get cell size on gather_lev - const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(gather_lev,0)); // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. @@ -2783,7 +2780,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, const amrex::ParticleReal Bz_external_particle = m_B_external_particle[2]; // Lower corner of tile box physical domain (take into account Galilean shift) - const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt); const Dim3 lo = lbound(box); @@ -2791,9 +2788,6 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); amrex::Array4 const& ez_arr = ezfab->array(); @@ -2908,7 +2902,7 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, + dinv, xyzmin, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); } @@ -3008,7 +3002,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, if (np_to_push == 0) { return; } // Get cell size on gather_lev - const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(gather_lev,0)); // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. @@ -3035,7 +3029,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, const amrex::ParticleReal Bz_external_particle = m_B_external_particle[2]; // Lower corner of tile box physical domain (take into account Galilean shift) - const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, gather_lev, 0._rt); const Dim3 lo = lbound(box); @@ -3043,9 +3037,6 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; - amrex::Array4 const& ex_arr = exfab->array(); amrex::Array4 const& ey_arr = eyfab->array(); amrex::Array4 const& ez_arr = ezfab->array(); @@ -3157,9 +3148,9 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, amrex::ParticleReal dxp, dxp_save; amrex::ParticleReal dyp, dyp_save; amrex::ParticleReal dzp, dzp_save; - auto idxg2 = static_cast(1._rt/(dx[0]*dx[0])); - auto idyg2 = static_cast(1._rt/(dx[1]*dx[1])); - auto idzg2 = static_cast(1._rt/(dx[2]*dx[2])); + auto idxg2 = static_cast(dinv.x*dinv.x); + auto idyg2 = static_cast(dinv.y*dinv.y); + auto idzg2 = static_cast(dinv.z*dinv.z); amrex::ParticleReal step_norm = 1._prt; for (int iter=0; iter& dx = WarpX::CellSize(std::max(lev,0)); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(lev,0)); #ifdef AMREX_USE_OMP #pragma omp parallel @@ -369,7 +369,6 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const amrex::ParticleReal By_external_particle = m_B_external_particle[1]; const amrex::ParticleReal Bz_external_particle = m_B_external_particle[2]; - const std::array& xyzmin = WarpX::LowerCorner(box, lev, 0._rt); const Dim3 lo = lbound(box); @@ -377,8 +376,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; - const amrex::GpuArray dx_arr = {dx[0], dx[1], dx[2]}; - const amrex::GpuArray xyzmin_arr = {xyzmin[0], xyzmin[1], xyzmin[2]}; + const amrex::XDim3 xyzmin = WarpX::LowerCorner(box, lev, 0._rt); amrex::Array4 const& ex_arr = exfab.array(); amrex::Array4 const& ey_arr = eyfab.array(); @@ -445,7 +443,7 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, doGatherShapeN(xp, yp, zp, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, ex_type, ey_type, ez_type, bx_type, by_type, bz_type, - dx_arr, xyzmin_arr, lo, n_rz_azimuthal_modes, + dinv, xyzmin, lo, n_rz_azimuthal_modes, nox, galerkin_interpolation); [[maybe_unused]] const auto& getExternalEB_tmp = getExternalEB; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 4800d9e209f..bdce18b7b2b 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -396,7 +396,8 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, amrex::numParticlesOutOfRange(pti, range) == 0, "Particles shape does not fit within tile (CPU) or guard cells (GPU) used for current deposition"); - const std::array& dx = WarpX::CellSize(std::max(depos_lev,0)); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(depos_lev,0)); + const amrex::ParticleReal q = this->charge; WARPX_PROFILE_VAR_NS("WarpXParticleContainer::DepositCurrent::Sorting", blp_sort); @@ -467,7 +468,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, // Note that this includes guard cells since it is after tilebox.ngrow const Dim3 lo = lbound(tilebox); // Take into account Galilean shift - const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev, 0.5_rt*dt); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, depos_lev, 0.5_rt*dt); if (WarpX::current_deposition_algo == CurrentDepositionAlgo::Esirkepov || WarpX::current_deposition_algo == CurrentDepositionAlgo::Villasenor) { @@ -544,28 +545,28 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, doDepositionSharedShapeN<1>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size); } else if (WarpX::nox == 2){ doDepositionSharedShapeN<2>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size); } else if (WarpX::nox == 3){ doDepositionSharedShapeN<3>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size); } else if (WarpX::nox == 4){ doDepositionSharedShapeN<4>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size); } @@ -580,25 +581,25 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, doEsirkepovDepositionShapeN<1>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doEsirkepovDepositionShapeN<2>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doEsirkepovDepositionShapeN<3>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 4){ doEsirkepovDepositionShapeN<4>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } } else if (push_type == PushType::Implicit) { @@ -625,7 +626,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doChargeConservingDepositionShapeNImplicit<2>( @@ -633,7 +634,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doChargeConservingDepositionShapeNImplicit<3>( @@ -641,7 +642,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 4){ doChargeConservingDepositionShapeNImplicit<4>( @@ -649,7 +650,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } } @@ -678,7 +679,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doVillasenorDepositionShapeNImplicit<2>( @@ -686,7 +687,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doVillasenorDepositionShapeNImplicit<3>( @@ -694,7 +695,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 4){ doVillasenorDepositionShapeNImplicit<4>( @@ -702,7 +703,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, GetPosition, wp.dataPtr() + offset, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dx, xyzmin, lo, q, + jx_arr, jy_arr, jz_arr, np_to_deposit, dt, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } } @@ -717,25 +718,25 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, doVayDepositionShapeN<1>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doVayDepositionShapeN<2>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doVayDepositionShapeN<3>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 4){ doVayDepositionShapeN<4>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dx, xyzmin, lo, q, + jx_fab, jy_fab, jz_fab, np_to_deposit, dt, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } } else { // Direct deposition @@ -744,25 +745,25 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, doDepositionShapeN<1>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doDepositionShapeN<2>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doDepositionShapeN<3>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 4){ doDepositionShapeN<4>( GetPosition, wp.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, relative_time, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } } else if (push_type == PushType::Implicit) { @@ -775,7 +776,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 2){ doDepositionShapeNImplicit<2>( @@ -783,7 +784,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 3){ doDepositionShapeNImplicit<3>( @@ -791,7 +792,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } else if (WarpX::nox == 4){ doDepositionShapeNImplicit<4>( @@ -799,7 +800,7 @@ WarpXParticleContainer::DepositCurrent (WarpXParIter& pti, uxp_n.dataPtr() + offset, uyp_n.dataPtr() + offset, uzp_n.dataPtr() + offset, uxp.dataPtr() + offset, uyp.dataPtr() + offset, uzp.dataPtr() + offset, ion_lev, - jx_fab, jy_fab, jz_fab, np_to_deposit, dx, + jx_fab, jy_fab, jz_fab, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes); } } @@ -931,7 +932,6 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, "Particles shape does not fit within tile (CPU) or guard cells (GPU) used for charge deposition"); amrex::ignore_unused(range); // In case the assertion isn't compiled - const std::array& dx = WarpX::CellSize(std::max(depos_lev,0)); const Real q = this->charge; WARPX_PROFILE_VAR_NS("WarpXParticleContainer::DepositCharge::Sorting", blp_sort); @@ -981,12 +981,13 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, // Take into account Galilean shift const Real dt = warpx.getdt(lev); const amrex::Real time_shift_delta = (icomp == 0 ? 0.0_rt : dt); - const std::array& xyzmin = WarpX::LowerCorner( - tilebox, depos_lev, time_shift_delta); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, depos_lev, time_shift_delta); // Indices of the lower bound const Dim3 lo = lbound(tilebox); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(depos_lev,0)); + // HACK - sort particles by bin here. WARPX_PROFILE_VAR_START(blp_sort); amrex::DenseBins bins; @@ -1088,25 +1089,25 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, if (WarpX::nox == 1){ doChargeDepositionSharedShapeN<1>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, ix_type, np_to_deposit, dx, xyzmin, lo, q, + rho_fab, ix_type, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size, WarpX::shared_tilesize); } else if (WarpX::nox == 2){ doChargeDepositionSharedShapeN<2>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, ix_type, np_to_deposit, dx, xyzmin, lo, q, + rho_fab, ix_type, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size, WarpX::shared_tilesize); } else if (WarpX::nox == 3){ doChargeDepositionSharedShapeN<3>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, ix_type, np_to_deposit, dx, xyzmin, lo, q, + rho_fab, ix_type, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size, WarpX::shared_tilesize); } else if (WarpX::nox == 4){ doChargeDepositionSharedShapeN<4>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, ix_type, np_to_deposit, dx, xyzmin, lo, q, + rho_fab, ix_type, np_to_deposit, dinv, xyzmin, lo, q, WarpX::n_rz_azimuthal_modes, bins, box, geom, max_tbox_size, WarpX::shared_tilesize); @@ -1125,7 +1126,6 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, // note: this is smaller than rho->nGrowVect() for PSATD const amrex::IntVect& ng_rho = warpx.get_ng_depos_rho(); - const std::array& dx = WarpX::CellSize(std::max(depos_lev,0)); amrex::IntVect ref_ratio; if (lev == depos_lev) { ref_ratio = IntVect(AMREX_D_DECL(1, 1, 1 )); @@ -1150,7 +1150,8 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, // Take into account Galilean shift const amrex::Real dt = warpx.getdt(lev); const amrex::Real time_shift_delta = (icomp == 0 ? 0.0_rt : dt); - const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev, time_shift_delta); + const amrex::XDim3 xyzmin = WarpX::LowerCorner(tilebox, depos_lev, time_shift_delta); + const amrex::XDim3 dinv = WarpX::InvCellSize(std::max(depos_lev,0)); AMREX_ALWAYS_ASSERT(WarpX::nox == WarpX::noy); AMREX_ALWAYS_ASSERT(WarpX::nox == WarpX::noz); @@ -1158,7 +1159,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector const& wp, ablastr::particles::deposit_charge( pti, wp, this->charge, ion_lev, rho, local_rho[thread_num], - WarpX::noz, dx, xyzmin, WarpX::n_rz_azimuthal_modes, + WarpX::noz, dinv, xyzmin, WarpX::n_rz_azimuthal_modes, ng_rho, depos_lev, ref_ratio, offset, np_to_deposit, icomp, nc); diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index d2691685c53..2ef4ee55d6e 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -232,8 +232,8 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) for(amrex::MFIter mfi(mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi){ const amrex::Box& bx = mfi.tilebox(); // Get box lower and upper physical z bound, and dz - const amrex::Real zmin_box = WarpX::LowerCorner(bx, lev, 0._rt)[2]; - const amrex::Real zmax_box = WarpX::UpperCorner(bx, lev, 0._rt)[2]; + const amrex::Real zmin_box = WarpX::LowerCorner(bx, lev, 0._rt).z; + const amrex::Real zmax_box = WarpX::UpperCorner(bx, lev, 0._rt).z; const amrex::Real dz = WarpX::CellSize(lev)[2]; // Get box lower index in the z direction #if defined(WARPX_DIM_3D) diff --git a/Source/WarpX.H b/Source/WarpX.H index a0a1379d9e2..768f2d486dc 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -895,6 +895,7 @@ public: amrex::Vector& output_geom ) const; static std::array CellSize (int lev); + static amrex::XDim3 InvCellSize (int lev); static amrex::RealBox getRealBox(const amrex::Box& bx, int lev); /** @@ -905,7 +906,7 @@ public: * (when v_galilean is not zero) * \return An array of the position coordinates */ - static std::array LowerCorner (const amrex::Box& bx, int lev, amrex::Real time_shift_delta); + static amrex::XDim3 LowerCorner (const amrex::Box& bx, int lev, amrex::Real time_shift_delta); /** * \brief Return the upper corner of the box in real units. * \param bx The box @@ -914,7 +915,7 @@ public: * (when v_galilean is not zero) * \return An array of the position coordinates */ - static std::array UpperCorner (const amrex::Box& bx, int lev, amrex::Real time_shift_delta); + static amrex::XDim3 UpperCorner (const amrex::Box& bx, int lev, amrex::Real time_shift_delta); static amrex::IntVect RefRatio (int lev); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index ef81aef4482..d1f9eb16f47 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -2846,6 +2846,13 @@ WarpX::CellSize (int lev) #endif } +amrex::XDim3 +WarpX::InvCellSize (int lev) +{ + std::array dx = WarpX::CellSize(lev); + return {1._rt/dx[0], 1._rt/dx[1], 1._rt/dx[2]}; +} + amrex::RealBox WarpX::getRealBox(const Box& bx, int lev) { @@ -2854,13 +2861,13 @@ WarpX::getRealBox(const Box& bx, int lev) return( grid_box ); } -std::array +amrex::XDim3 WarpX::LowerCorner(const Box& bx, const int lev, const amrex::Real time_shift_delta) { auto & warpx = GetInstance(); const RealBox grid_box = getRealBox( bx, lev ); - const Real* xyzmin = grid_box.lo(); + const Real* grid_min = grid_box.lo(); const amrex::Real cur_time = warpx.gett_new(lev); const amrex::Real time_shift = (cur_time + time_shift_delta - warpx.time_of_last_gal_shift); @@ -2869,23 +2876,23 @@ WarpX::LowerCorner(const Box& bx, const int lev, const amrex::Real time_shift_de warpx.m_v_galilean[2]*time_shift }; #if defined(WARPX_DIM_3D) - return { xyzmin[0] + galilean_shift[0], xyzmin[1] + galilean_shift[1], xyzmin[2] + galilean_shift[2] }; + return { grid_min[0] + galilean_shift[0], grid_min[1] + galilean_shift[1], grid_min[2] + galilean_shift[2] }; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - return { xyzmin[0] + galilean_shift[0], std::numeric_limits::lowest(), xyzmin[1] + galilean_shift[2] }; + return { grid_min[0] + galilean_shift[0], std::numeric_limits::lowest(), grid_min[1] + galilean_shift[2] }; #elif defined(WARPX_DIM_1D_Z) - return { std::numeric_limits::lowest(), std::numeric_limits::lowest(), xyzmin[0] + galilean_shift[2] }; + return { std::numeric_limits::lowest(), std::numeric_limits::lowest(), grid_min[0] + galilean_shift[2] }; #endif } -std::array +amrex::XDim3 WarpX::UpperCorner(const Box& bx, const int lev, const amrex::Real time_shift_delta) { auto & warpx = GetInstance(); const RealBox grid_box = getRealBox( bx, lev ); - const Real* xyzmax = grid_box.hi(); + const Real* grid_max = grid_box.hi(); const amrex::Real cur_time = warpx.gett_new(lev); const amrex::Real time_shift = (cur_time + time_shift_delta - warpx.time_of_last_gal_shift); @@ -2894,13 +2901,13 @@ WarpX::UpperCorner(const Box& bx, const int lev, const amrex::Real time_shift_de warpx.m_v_galilean[2]*time_shift }; #if defined(WARPX_DIM_3D) - return { xyzmax[0] + galilean_shift[0], xyzmax[1] + galilean_shift[1], xyzmax[2] + galilean_shift[2] }; + return { grid_max[0] + galilean_shift[0], grid_max[1] + galilean_shift[1], grid_max[2] + galilean_shift[2] }; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - return { xyzmax[0] + galilean_shift[0], std::numeric_limits::max(), xyzmax[1] + galilean_shift[1] }; + return { grid_max[0] + galilean_shift[0], std::numeric_limits::max(), grid_max[1] + galilean_shift[1] }; #elif defined(WARPX_DIM_1D_Z) - return { std::numeric_limits::max(), std::numeric_limits::max(), xyzmax[0] + galilean_shift[0] }; + return { std::numeric_limits::max(), std::numeric_limits::max(), grid_max[0] + galilean_shift[0] }; #endif } diff --git a/Source/ablastr/particles/DepositCharge.H b/Source/ablastr/particles/DepositCharge.H index ff3741a7a43..75e3bca170d 100644 --- a/Source/ablastr/particles/DepositCharge.H +++ b/Source/ablastr/particles/DepositCharge.H @@ -34,7 +34,7 @@ namespace ablastr::particles * \param rho MultiFab of the charge density * \param local_rho temporary FArrayBox for deposition with OpenMP * \param particle_shape shape factor in each direction - * \param dx cell spacing at level lev + * \param dinv cell spacing inverses at level lev * \param xyzmin lo corner of the current tile in physical coordinates. * \param n_rz_azimuthal_modes number of azimuthal modes in use, irrelevant outside RZ geometry (default: 0) * \param num_rho_deposition_guards number of ghost cells to use for rho (default: rho.nGrowVect()) @@ -54,8 +54,8 @@ deposit_charge (typename T_PC::ParIterType& pti, amrex::MultiFab* rho, amrex::FArrayBox& local_rho, int const particle_shape, - std::array const & dx, - std::array const & xyzmin, + const amrex::XDim3 dinv, + const amrex::XDim3 xyzmin, int const n_rz_azimuthal_modes = 0, std::optional num_rho_deposition_guards = std::nullopt, std::optional depos_lev = std::nullopt, @@ -175,19 +175,19 @@ deposit_charge (typename T_PC::ParIterType& pti, if (nox == 1){ doChargeDepositionShapeN<1>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge, + rho_fab, np_to_deposit.value(), dinv, xyzmin, lo, charge, n_rz_azimuthal_modes); } else if (nox == 2){ doChargeDepositionShapeN<2>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge, + rho_fab, np_to_deposit.value(), dinv, xyzmin, lo, charge, n_rz_azimuthal_modes); } else if (nox == 3){ doChargeDepositionShapeN<3>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge, + rho_fab, np_to_deposit.value(), dinv, xyzmin, lo, charge, n_rz_azimuthal_modes); } else if (nox == 4){ doChargeDepositionShapeN<4>(GetPosition, wp.dataPtr()+offset, ion_lev, - rho_fab, np_to_deposit.value(), dx, xyzmin, lo, charge, + rho_fab, np_to_deposit.value(), dinv, xyzmin, lo, charge, n_rz_azimuthal_modes); } ABLASTR_PROFILE_VAR_STOP(blp_ppc_chd);