diff options
Diffstat (limited to 'Source')
15 files changed, 254 insertions, 241 deletions
diff --git a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp index c0ecd5934..0b50e2d79 100644 --- a/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp +++ b/Source/Diagnostics/ComputeDiagFunctors/ParticleReductionFunctor.cpp @@ -92,7 +92,7 @@ ParticleReductionFunctor::operator() (amrex::MultiFab& mf_dst, const int dcomp, amrex::Real value; if ((do_filter) && (!filter_fn(xw, yw, zw, ux, uy, uz))) value = 0._rt; else value = map_fn(xw, yw, zw, ux, uy, uz); - amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), p.rdata(PIdx::w) * value); + amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), (amrex::Real)(p.rdata(PIdx::w) * value)); }); if (m_do_average) { amrex::MultiFab ppc_mf(warpx.boxArray(m_lev), warpx.DistributionMap(m_lev), 1, ng); @@ -134,7 +134,7 @@ ParticleReductionFunctor::operator() (amrex::MultiFab& mf_dst, const int dcomp, amrex::Real filter; if ((do_filter) && (!filter_fn(xw, yw, zw, ux, uy, uz))) filter = 0._rt; else filter = 1._rt; - amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), p.rdata(PIdx::w) * filter); + amrex::Gpu::Atomic::AddNoRet(&out_array(ii, jj, kk, 0), (amrex::Real)(p.rdata(PIdx::w) * filter)); }); // Divide value by number of particles for average. Set average to zero if there are no particles for (amrex::MFIter mfi(red_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) diff --git a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp index 65f539553..c63071faa 100644 --- a/Source/Diagnostics/ReducedDiags/FieldProbe.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldProbe.cpp @@ -545,7 +545,7 @@ void FieldProbe::ComputeDiags (int step) temp_interp_order, false); //Calculate the Poynting Vector S - amrex::Real const sraw[3]{ + amrex::ParticleReal const sraw[3]{ Exp * Bzp - Ezp * Byp, Ezp * Bxp - Exp * Bzp, Exp * Byp - Eyp * Bxp diff --git a/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp b/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp index 415083afd..f54f9bf37 100644 --- a/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp +++ b/Source/Diagnostics/ReducedDiags/ParticleHistogram.cpp @@ -198,10 +198,10 @@ void ParticleHistogram::ComputeDiags (int step) auto const GetPosition = GetParticlePosition(pti); auto & attribs = pti.GetAttribs(); - Real* const AMREX_RESTRICT d_w = attribs[PIdx::w].dataPtr(); - Real* const AMREX_RESTRICT d_ux = attribs[PIdx::ux].dataPtr(); - Real* const AMREX_RESTRICT d_uy = attribs[PIdx::uy].dataPtr(); - Real* const AMREX_RESTRICT d_uz = attribs[PIdx::uz].dataPtr(); + ParticleReal* const AMREX_RESTRICT d_w = attribs[PIdx::w].dataPtr(); + ParticleReal* const AMREX_RESTRICT d_ux = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT d_uy = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT d_uz = attribs[PIdx::uz].dataPtr(); long const np = pti.numParticles(); @@ -211,7 +211,7 @@ void ParticleHistogram::ComputeDiags (int step) { amrex::ParticleReal x, y, z; GetPosition(i, x, y, z); - auto const w = d_w[i]; + auto const w = (amrex::Real)d_w[i]; auto const ux = d_ux[i] / PhysConst::c; auto const uy = d_uy[i] / PhysConst::c; auto const uz = d_uz[i] / PhysConst::c; @@ -222,7 +222,6 @@ void ParticleHistogram::ComputeDiags (int step) return; // continue function if particle is not filtered out auto const f = fun_partparser(t, x, y, z, ux, uy, uz); - // determine particle bin int const bin = int(Math::floor((f-bin_min)/bin_size)); if ( bin<0 || bin>=num_bins ) return; // discard if out-of-range diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index b980e9403..41e11ea63 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -119,6 +119,20 @@ WarpX::PrintMainPICparameters () amrex::Print() << "--------------------------- MAIN EM PIC PARAMETERS ----------------------------\n"; amrex::Print() << "-------------------------------------------------------------------------------\n"; + // print warpx build information + if constexpr (std::is_same<Real, float>::value) { + amrex::Print() << "Precision: | SINGLE" << "\n"; + } + else { + amrex::Print() << "Precision: | DOUBLE" << "\n"; + } + if constexpr (std::is_same<ParticleReal, float>::value) { + amrex::Print() << "Particle precision: | SINGLE" << "\n"; + } + else { + amrex::Print() << "Particle precision: | DOUBLE" << "\n"; + } + // Print geometry dimensionality amrex::ParmParse pp_geometry("geometry"); std::string dims; diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index 658218962..3387f606a 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -74,9 +74,9 @@ public: index_type const* AMREX_RESTRICT I2, SoaData_type soa_1, SoaData_type soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, - amrex::Real const q1, amrex::Real const q2, - amrex::Real const m1, amrex::Real const m2, - amrex::Real const dt, amrex::Real const dV, + amrex::ParticleReal const q1, amrex::ParticleReal const q2, + amrex::ParticleReal const m1, amrex::ParticleReal const m2, + amrex::ParticleReal const dt, amrex::ParticleReal const dV, index_type const /*cell_start_pair*/, index_type* /*p_mask*/, index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/, amrex::ParticleReal* /*p_pair_reaction_weight*/, @@ -85,12 +85,12 @@ public: ElasticCollisionPerez( I1s, I1e, I2s, I2e, I1, I2, soa_1, soa_2, - q1, q2, m1, m2, amrex::Real(-1.0), amrex::Real(-1.0), + q1, q2, m1, m2, -1.0_prt, -1.0_prt, dt, m_CoulombLog, dV, engine ); } private: - amrex::Real m_CoulombLog; + amrex::ParticleReal m_CoulombLog; }; #endif // PAIRWISE_COULOMB_COLLISION_FUNC_H_ diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H index 38ddc5481..fb7dc4da1 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H @@ -106,7 +106,8 @@ void UpdateMomentumPerezElastic ( ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_R(1.0) ); // Compute the minimal impact parameter - T_R bmin = amrex::max(PhysConst::hbar*MathConst::pi/p1sm,b0); + constexpr T_R hbar_pi = static_cast<T_R>(PhysConst::hbar*MathConst::pi); + const T_R bmin = amrex::max(hbar_pi/p1sm, b0); // Compute the Coulomb log lnLmd lnLmd = amrex::max( T_R(2.0), diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H index c59ad1a12..30c054bc4 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H @@ -58,10 +58,10 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part const amrex::Real& dt, const amrex::Real& dV, const int& pair_index, index_type* AMREX_RESTRICT p_mask, amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, - const amrex::Real& fusion_multiplier, + const amrex::ParticleReal& fusion_multiplier, const int& multiplier_ratio, - const amrex::Real& probability_threshold, - const amrex::Real& probability_target_value, + const amrex::ParticleReal& probability_threshold, + const amrex::ParticleReal& probability_target_value, const NuclearFusionType& fusion_type, const amrex::RandomEngine& engine) { diff --git a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H index cf04d86fa..456fdbb5d 100644 --- a/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H +++ b/Source/Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H @@ -175,7 +175,7 @@ public: return 0; const auto is_out = pxr_bw::evolve_optical_depth< - amrex::Real, + amrex::ParticleReal, BW_dndt_table_view, pxr_p::unit_system::SI>( energy, chi_phot, dt, opt_depth, m_table_view); diff --git a/Source/Particles/Filter/FilterFunctors.H b/Source/Particles/Filter/FilterFunctors.H index c752ea272..00f960ef9 100644 --- a/Source/Particles/Filter/FilterFunctors.H +++ b/Source/Particles/Filter/FilterFunctors.H @@ -125,7 +125,7 @@ struct ParserFilter uz /= m_mass; } // ux, uy, uz are now in beta*gamma - if ( m_function_partparser(m_t,x,y,z,ux,uy,uz) ) return 1; + if ( m_function_partparser(m_t,x,y,z,ux,uy,uz)) return 1; else return 0; } private: @@ -135,9 +135,9 @@ public: /** Parser function with 7 input variables, t,x,y,z,ux,uy,uz */ amrex::ParserExecutor<7> const m_function_partparser; /** Store physical time. */ - amrex::Real m_t; + amrex::ParticleReal m_t; /** Mass of particle species */ - amrex::Real m_mass; + amrex::ParticleReal m_mass; /** keep track of momentum units particles will come in with **/ InputUnits m_units; }; diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H index 3c1dfdf60..954b94fc7 100644 --- a/Source/Particles/Gather/GetExternalFields.H +++ b/Source/Particles/Gather/GetExternalFields.H @@ -91,9 +91,9 @@ struct GetExternalEBField lab_time = m_gamma_boost*m_time + m_uz_boost*z*inv_c2; z = m_gamma_boost*z + m_uz_boost*m_time; } - Ex = m_Exfield_partparser(x, y, z, lab_time); - Ey = m_Eyfield_partparser(x, y, z, lab_time); - Ez = m_Ezfield_partparser(x, y, z, lab_time); + Ex = m_Exfield_partparser((amrex::Real) x, (amrex::Real) y, (amrex::Real) z, lab_time); + Ey = m_Eyfield_partparser((amrex::Real) x, (amrex::Real) y, (amrex::Real) z, lab_time); + Ez = m_Ezfield_partparser((amrex::Real) x, (amrex::Real) y, (amrex::Real) z, lab_time); } if (m_Btype == Constant) diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index 12d7e78cd..fbe2b9710 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -453,7 +453,7 @@ LaserParticleContainer::InitData (int lev) BoxArray plane_ba { Box {IntVect(0), IntVect(0)} }; #endif - amrex::Vector<amrex::Real> particle_x, particle_y, particle_z, particle_w; + amrex::Vector<amrex::ParticleReal> particle_x, particle_y, particle_z, particle_w; const DistributionMapping plane_dm {plane_ba, nprocs}; const Vector<int>& procmap = plane_dm.ProcessorMap(); @@ -506,9 +506,9 @@ LaserParticleContainer::InitData (int lev) } } const int np = particle_z.size(); - amrex::Vector<amrex::Real> particle_ux(np, 0.0); - amrex::Vector<amrex::Real> particle_uy(np, 0.0); - amrex::Vector<amrex::Real> particle_uz(np, 0.0); + amrex::Vector<amrex::ParticleReal> particle_ux(np, 0.0); + amrex::Vector<amrex::ParticleReal> particle_uy(np, 0.0); + amrex::Vector<amrex::ParticleReal> particle_uz(np, 0.0); if (Verbose()) amrex::Print() << Utils::TextMsg::Info("Adding laser particles"); // Add particles on level 0. They will be redistributed afterwards diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index bc0366adc..47ece0d31 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -193,8 +193,8 @@ public: */ void AddPlasmaFlux (amrex::Real dt); - void MapParticletoBoostedFrame (amrex::Real& x, amrex::Real& y, amrex::Real& z, - amrex::Real& ux, amrex::Real& uy, amrex::Real& uz); + void MapParticletoBoostedFrame (amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, + amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz); void AddGaussianBeam ( const amrex::Real x_m, const amrex::Real y_m, const amrex::Real z_m, @@ -210,9 +210,9 @@ public: amrex::ParticleReal z_shift); void CheckAndAddParticle ( - amrex::Real x, amrex::Real y, amrex::Real z, - amrex::Real ux, amrex::Real uy, amrex::Real uz, - amrex::Real weight, + amrex::ParticleReal x, amrex::ParticleReal y, amrex::ParticleReal z, + amrex::ParticleReal ux, amrex::ParticleReal uy, amrex::ParticleReal uz, + amrex::ParticleReal weight, amrex::Gpu::HostVector<amrex::ParticleReal>& particle_x, amrex::Gpu::HostVector<amrex::ParticleReal>& particle_y, amrex::Gpu::HostVector<amrex::ParticleReal>& particle_z, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ccd5edf6c..20632ce1f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -137,6 +137,20 @@ namespace return z0; } + struct PDim3 { + ParticleReal x, y, z; + + AMREX_GPU_HOST_DEVICE + PDim3(const PDim3&) = default; + AMREX_GPU_HOST_DEVICE + PDim3(const amrex::XDim3& a) : x(a.x), y(a.y), z(a.z) {} + + AMREX_GPU_HOST_DEVICE + PDim3& operator=(const PDim3&) = default; + AMREX_GPU_HOST_DEVICE + PDim3& operator=(const amrex::XDim3 &a) { x = a.x; y = a.y; z = a.z; return *this; } + }; + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE XDim3 getCellCoords (const GpuArray<Real, AMREX_SPACEDIM>& lo_corner, const GpuArray<Real, AMREX_SPACEDIM>& dx, @@ -377,7 +391,7 @@ void PhysicalParticleContainer::InitData () } void PhysicalParticleContainer::MapParticletoBoostedFrame ( - Real& x, Real& y, Real& z, Real& ux, Real& uy, Real& uz) + ParticleReal& x, ParticleReal& y, ParticleReal& z, ParticleReal& ux, ParticleReal& uy, ParticleReal& uz) { // Map the particles from the lab frame to the boosted frame. // This boosts the particle to the lab frame and calculates @@ -386,28 +400,28 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame ( // For now, start with the assumption that this will only happen // at the start of the simulation. - const Real t_lab = 0._rt; + const ParticleReal t_lab = 0._prt; - const Real uz_boost = WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; + const ParticleReal uz_boost = WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; // tpr is the particle's time in the boosted frame - Real tpr = WarpX::gamma_boost*t_lab - uz_boost*z/(PhysConst::c*PhysConst::c); + ParticleReal tpr = WarpX::gamma_boost*t_lab - uz_boost*z/(PhysConst::c*PhysConst::c); // The particle's transformed location in the boosted frame - Real xpr = x; - Real ypr = y; - Real zpr = WarpX::gamma_boost*z - uz_boost*t_lab; + ParticleReal xpr = x; + ParticleReal ypr = y; + ParticleReal zpr = WarpX::gamma_boost*z - uz_boost*t_lab; // transform u and gamma to the boosted frame - Real gamma_lab = std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)/(PhysConst::c*PhysConst::c)); + ParticleReal gamma_lab = std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)/(PhysConst::c*PhysConst::c)); // ux = ux; // uy = uy; uz = WarpX::gamma_boost*uz - uz_boost*gamma_lab; - Real gammapr = std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)/(PhysConst::c*PhysConst::c)); + ParticleReal gammapr = std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)/(PhysConst::c*PhysConst::c)); - Real vxpr = ux/gammapr; - Real vypr = uy/gammapr; - Real vzpr = uz/gammapr; + ParticleReal vxpr = ux/gammapr; + ParticleReal vypr = uy/gammapr; + ParticleReal vzpr = uz/gammapr; if (do_backward_propagation){ uz = -uz; @@ -622,9 +636,9 @@ PhysicalParticleContainer::AddPlasmaFromFile(ParticleReal q_tot, void PhysicalParticleContainer::CheckAndAddParticle ( - Real x, Real y, Real z, - Real ux, Real uy, Real uz, - Real weight, + ParticleReal x, ParticleReal y, ParticleReal z, + ParticleReal ux, ParticleReal uy, ParticleReal uz, + ParticleReal weight, Gpu::HostVector<ParticleReal>& particle_x, Gpu::HostVector<ParticleReal>& particle_y, Gpu::HostVector<ParticleReal>& particle_z, @@ -1496,32 +1510,34 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) const XDim3 r = inj_pos->getPositionUnitBox(i_part, lrrfac, engine); auto pos = getCellCoords(overlap_corner, dx, r, iv); + auto ppos = PDim3(pos); // inj_mom would typically be InjectorMomentumGaussianFlux XDim3 u; u = inj_mom->getMomentum(pos.x, pos.y, pos.z, engine); + auto pu = PDim3(u); - u.x *= PhysConst::c; - u.y *= PhysConst::c; - u.z *= PhysConst::c; + pu.x *= PhysConst::c; + pu.y *= PhysConst::c; + pu.z *= PhysConst::c; const amrex::Real t_fract = amrex::Random(engine)*dt; - UpdatePosition(pos.x, pos.y, pos.z, u.x, u.y, u.z, t_fract); + UpdatePosition(ppos.x, ppos.y, ppos.z, pu.x, pu.y, pu.z, t_fract); #if defined(WARPX_DIM_3D) - if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) { + if (!tile_realbox.contains(XDim3{ppos.x,ppos.y,ppos.z})) { p.id() = -1; continue; } #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(k); - if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { + if (!tile_realbox.contains(XDim3{ppos.x,ppos.z,0.0_prt})) { p.id() = -1; continue; } #else amrex::ignore_unused(j,k); - if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { + if (!tile_realbox.contains(XDim3{ppos.z,0.0_prt,0.0_prt})) { p.id() = -1; continue; } @@ -1529,8 +1545,8 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) // Save the x and y values to use in the insideBounds checks. // This is needed with WARPX_DIM_RZ since x and y are modified. - Real xb = pos.x; - Real yb = pos.y; + Real xb = ppos.x; + Real yb = ppos.y; #ifdef WARPX_DIM_RZ // Replace the x and y, setting an angle theta. @@ -1539,25 +1555,25 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) if (nmodes == 1) { // With only 1 mode, the angle doesn't matter so // choose it randomly. - theta = 2._rt*MathConst::pi*amrex::Random(engine); + theta = 2._prt*MathConst::pi*amrex::Random(engine); } else { - theta = 2._rt*MathConst::pi*r.y; + theta = 2._prt*MathConst::pi*r.y; } Real const cos_theta = std::cos(theta); Real const sin_theta = std::sin(theta); // Rotate the position - pos.x = xb*cos_theta; - pos.y = xb*sin_theta; + ppos.x = xb*cos_theta; + ppos.y = xb*sin_theta; if (loc_flux_normal_axis != 2) { // Rotate the momentum // This because, when the flux direction is e.g. "r" // the `inj_mom` objects generates a v*Gaussian distribution // along the Cartesian "x" directionm by default. This // needs to be rotated along "r". - Real ur = u.x; - Real ut = u.y; - u.x = cos_theta*ur - sin_theta*ut; - u.y = sin_theta*ur + cos_theta*ut; + Real ur = pu.x; + Real ut = pu.y; + pu.x = cos_theta*ur - sin_theta*ut; + pu.y = sin_theta*ur + cos_theta*ut; } #endif @@ -1566,12 +1582,12 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) // xmin, xmax, ymin, ymax, zmin, zmax, go to // the next generated particle. - if (!inj_pos->insideBounds(xb, yb, pos.z)) { + if (!inj_pos->insideBounds(xb, yb, ppos.z)) { p.id() = -1; continue; } - Real dens = inj_rho->getDensity(pos.x, pos.y, pos.z); + Real dens = inj_rho->getDensity(ppos.x, ppos.y, ppos.z); // Remove particle if density below threshold if ( dens < density_min ){ @@ -1609,28 +1625,28 @@ PhysicalParticleContainer::AddPlasmaFlux (amrex::Real dt) } else { // This is not correct since it might shift the particle // out of the local grid - pos.x = std::sqrt(xb*rmax); + ppos.x = std::sqrt(xb*rmax); weight *= dx[0]; } } #endif pa[PIdx::w ][ip] = weight; - pa[PIdx::ux][ip] = u.x; - pa[PIdx::uy][ip] = u.y; - pa[PIdx::uz][ip] = u.z; + pa[PIdx::ux][ip] = pu.x; + pa[PIdx::uy][ip] = pu.y; + pa[PIdx::uz][ip] = pu.z; #if defined(WARPX_DIM_3D) - p.pos(0) = pos.x; - p.pos(1) = pos.y; - p.pos(2) = pos.z; + p.pos(0) = ppos.x; + p.pos(1) = ppos.y; + p.pos(2) = ppos.z; #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) #ifdef WARPX_DIM_RZ pa[PIdx::theta][ip] = theta; #endif p.pos(0) = xb; - p.pos(1) = pos.z; + p.pos(1) = ppos.z; #else - p.pos(0) = pos.z; + p.pos(0) = ppos.z; #endif } }); @@ -2363,22 +2379,22 @@ PhysicalParticleContainer::GetParticleSlice ( const auto GetPosition = GetParticlePosition(pti); auto& attribs = pti.GetAttribs(); - Real* const AMREX_RESTRICT wpnew = attribs[PIdx::w].dataPtr(); - Real* const AMREX_RESTRICT uxpnew = attribs[PIdx::ux].dataPtr(); - Real* const AMREX_RESTRICT uypnew = attribs[PIdx::uy].dataPtr(); - Real* const AMREX_RESTRICT uzpnew = attribs[PIdx::uz].dataPtr(); + ParticleReal* const AMREX_RESTRICT wpnew = attribs[PIdx::w].dataPtr(); + ParticleReal* const AMREX_RESTRICT uxpnew = attribs[PIdx::ux].dataPtr(); + ParticleReal* const AMREX_RESTRICT uypnew = attribs[PIdx::uy].dataPtr(); + ParticleReal* const AMREX_RESTRICT uzpnew = attribs[PIdx::uz].dataPtr(); - Real* const AMREX_RESTRICT + ParticleReal* const AMREX_RESTRICT xpold = tmp_particle_data[lev][index][TmpIdx::xold].dataPtr(); - Real* const AMREX_RESTRICT + ParticleReal* const AMREX_RESTRICT ypold = tmp_particle_data[lev][index][TmpIdx::yold].dataPtr(); - Real* const AMREX_RESTRICT + ParticleReal* const AMREX_RESTRICT zpold = tmp_particle_data[lev][index][TmpIdx::zold].dataPtr(); - Real* const AMREX_RESTRICT + ParticleReal* const AMREX_RESTRICT uxpold = tmp_particle_data[lev][index][TmpIdx::uxold].dataPtr(); - Real* const AMREX_RESTRICT + ParticleReal* const AMREX_RESTRICT uypold = tmp_particle_data[lev][index][TmpIdx::uyold].dataPtr(); - Real* const AMREX_RESTRICT + ParticleReal* const AMREX_RESTRICT uzpold = tmp_particle_data[lev][index][TmpIdx::uzold].dataPtr(); const long np = pti.numParticles(); @@ -2418,19 +2434,19 @@ PhysicalParticleContainer::GetParticleSlice ( amrex::Real betaboost = WarpX::beta_boost; amrex::Real Phys_c = PhysConst::c; - Real* const AMREX_RESTRICT diag_wp = + ParticleReal* const AMREX_RESTRICT diag_wp = diagnostic_particles[lev][index].GetRealData(DiagIdx::w).data(); - Real* const AMREX_RESTRICT diag_xp = + ParticleReal* const AMREX_RESTRICT diag_xp = diagnostic_particles[lev][index].GetRealData(DiagIdx::x).data(); - Real* const AMREX_RESTRICT diag_yp = + ParticleReal* const AMREX_RESTRICT diag_yp = diagnostic_particles[lev][index].GetRealData(DiagIdx::y).data(); - Real* const AMREX_RESTRICT diag_zp = + ParticleReal* const AMREX_RESTRICT diag_zp = diagnostic_particles[lev][index].GetRealData(DiagIdx::z).data(); - Real* const AMREX_RESTRICT diag_uxp = + ParticleReal* const AMREX_RESTRICT diag_uxp = diagnostic_particles[lev][index].GetRealData(DiagIdx::ux).data(); - Real* const AMREX_RESTRICT diag_uyp = + ParticleReal* const AMREX_RESTRICT diag_uyp = diagnostic_particles[lev][index].GetRealData(DiagIdx::uy).data(); - Real* const AMREX_RESTRICT diag_uzp = + ParticleReal* const AMREX_RESTRICT diag_uzp = diagnostic_particles[lev][index].GetRealData(DiagIdx::uz).data(); // Copy particle data to diagnostic particle array on the GPU diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index f5f434932..9360c14d4 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -247,8 +247,56 @@ amrex::ParserExecutor<N> compileParser (amrex::Parser const* parser) * \param[in] str name of the parameter to read * \param[out] val where the value queried and parsed is stored, either a scalar or vector */ -int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, float& val); -int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, double& val); +template <typename T> +int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, T& val) +{ + // call amrex::ParmParse::query, check if the user specified str. + std::string tmp_str; + int is_specified = a_pp.query(str, tmp_str); + if (is_specified) + { + // If so, create a parser object and apply it to the value provided by the user. + std::string str_val; + Store_parserString(a_pp, str, str_val); + + auto parser = makeParser(str_val, {}); + + if (std::is_same<T, int>::value) { + + val = safeCastToInt(std::round(parser.compileHost<0>()()), str); + } + else { + val = static_cast<T>(parser.compileHost<0>()()); + } + } + // return the same output as amrex::ParmParse::query + return is_specified; +} + +template <typename T> +int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<T>& val) +{ + // call amrex::ParmParse::query, check if the user specified str. + std::vector<std::string> tmp_str_arr; + int is_specified = a_pp.queryarr(str, tmp_str_arr); + if (is_specified) + { + // If so, create parser objects and apply them to the values provided by the user. + int const n = static_cast<int>(tmp_str_arr.size()); + val.resize(n); + for (int i=0 ; i < n ; i++) { + auto parser = makeParser(tmp_str_arr[i], {}); + if (std::is_same<T, int>::value) { + val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); + } + else { + val[i] = static_cast<T>(parser.compileHost<0>()()); + } + } + } + // return the same output as amrex::ParmParse::query + return is_specified; +} /** Similar to amrex::ParmParse::query, but also supports math expressions for the value. * @@ -265,13 +313,31 @@ int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, doubl * \param[in] num_val number of input values to use (optional with arrays, default is * amrex::ParmParse::LAST for reading until the last input value) */ -int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<amrex::Real>& val, - const int start_ix = amrex::ParmParse::FIRST, - const int num_val = amrex::ParmParse::LAST); -int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, int& val); -int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<int>& val, - const int start_ix = amrex::ParmParse::FIRST, - const int num_val = amrex::ParmParse::LAST); +template <typename T> +int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<T>& val, + const int start_ix, const int num_val) +{ + // call amrex::ParmParse::query, check if the user specified str. + std::vector<std::string> tmp_str_arr; + int is_specified = a_pp.queryarr(str, tmp_str_arr, start_ix, num_val); + if (is_specified) + { + // If so, create parser objects and apply them to the values provided by the user. + int const n = static_cast<int>(tmp_str_arr.size()); + val.resize(n); + for (int i=0 ; i < n ; i++) { + auto parser = makeParser(tmp_str_arr[i], {}); + if (std::is_same<T, int>::value) { + val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); + } + else { + val[i] = static_cast<T>(parser.compileHost<0>()()); + } + } + } + // return the same output as amrex::ParmParse::query + return is_specified; +} /** Similar to amrex::ParmParse::get, but also supports math expressions for the value. * @@ -284,8 +350,41 @@ int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, st * \param[in] str name of the parameter to read * \param[out] val where the value queried and parsed is stored */ -void getWithParser (const amrex::ParmParse& a_pp, char const * const str, float& val); -void getWithParser (const amrex::ParmParse& a_pp, char const * const str, double& val); +template <typename T> +void getWithParser (const amrex::ParmParse& a_pp, char const * const str, T& val) +{ + // If so, create a parser object and apply it to the value provided by the user. + std::string str_val; + Store_parserString(a_pp, str, str_val); + + auto parser = makeParser(str_val, {}); + if (std::is_same<T, int>::value) { + val = safeCastToInt(std::round(parser.compileHost<0>()()), str); + } + else { + val = static_cast<T>(parser.compileHost<0>()()); + } +} + +template <typename T> +void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<T>& val) +{ + // Create parser objects and apply them to the values provided by the user. + std::vector<std::string> tmp_str_arr; + a_pp.getarr(str, tmp_str_arr); + + int const n = static_cast<int>(tmp_str_arr.size()); + val.resize(n); + for (int i=0 ; i < n ; i++) { + auto parser = makeParser(tmp_str_arr[i], {}); + if (std::is_same<T, int>::value) { + val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); + } + else { + val[i] = static_cast<T>(parser.compileHost<0>()()); + } + } +} /** Similar to amrex::ParmParse::get, but also supports math expressions for the value. * @@ -302,13 +401,26 @@ void getWithParser (const amrex::ParmParse& a_pp, char const * const str, double * \param[in] num_val number of input values to use (optional with arrays, default is * amrex::ParmParse::LAST for reading until the last input value) */ -void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<amrex::Real>& val, - const int start_ix = amrex::ParmParse::FIRST, - const int num_val = amrex::ParmParse::LAST); -void getWithParser (const amrex::ParmParse& a_pp, char const * const str, int& val); -void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<int>& val, - const int start_ix = amrex::ParmParse::FIRST, - const int num_val = amrex::ParmParse::LAST); +template <typename T> +void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<T>& val, + const int start_ix, const int num_val) +{ + // Create parser objects and apply them to the values provided by the user. + std::vector<std::string> tmp_str_arr; + a_pp.getarr(str, tmp_str_arr, start_ix, num_val); + + int const n = static_cast<int>(tmp_str_arr.size()); + val.resize(n); + for (int i=0 ; i < n ; i++) { + auto parser = makeParser(tmp_str_arr[i], {}); + if (std::is_same<T, int>::value) { + val[i] = safeCastToInt(std::round(parser.compileHost<0>()()), str); + } + else { + val[i] = static_cast<T>(parser.compileHost<0>()()); + } + } +} namespace WarpXUtilStr { diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 9ab96873a..afe8b7daa 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -388,135 +388,6 @@ parseStringtoInt(std::string str, std::string name) return ival; } -// Overloads for float/double instead of amrex::Real to allow makeParser() to query for -// my_constants as double even in single precision mode -// Always parsing in double precision avoids potential overflows that may occur when parsing user's -// expressions because of the limited range of exponentials in single precision -int -queryWithParser (const amrex::ParmParse& a_pp, char const * const str, float& val) -{ - // call amrex::ParmParse::query, check if the user specified str. - std::string tmp_str; - int is_specified = a_pp.query(str, tmp_str); - if (is_specified) - { - // If so, create a parser object and apply it to the value provided by the user. - std::string str_val; - Store_parserString(a_pp, str, str_val); - val = static_cast<float>(parseStringtoReal(str_val)); - } - // return the same output as amrex::ParmParse::query - return is_specified; -} - -void -getWithParser (const amrex::ParmParse& a_pp, char const * const str, float& val) -{ - // If so, create a parser object and apply it to the value provided by the user. - std::string str_val; - Store_parserString(a_pp, str, str_val); - val = static_cast<float>(parseStringtoReal(str_val)); -} - -int -queryWithParser (const amrex::ParmParse& a_pp, char const * const str, double& val) -{ - // call amrex::ParmParse::query, check if the user specified str. - std::string tmp_str; - int is_specified = a_pp.query(str, tmp_str); - if (is_specified) - { - // If so, create a parser object and apply it to the value provided by the user. - std::string str_val; - Store_parserString(a_pp, str, str_val); - val = parseStringtoReal(str_val); - } - // return the same output as amrex::ParmParse::query - return is_specified; -} - -void -getWithParser (const amrex::ParmParse& a_pp, char const * const str, double& val) -{ - // If so, create a parser object and apply it to the value provided by the user. - std::string str_val; - Store_parserString(a_pp, str, str_val); - val = parseStringtoReal(str_val); -} - -int -queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<amrex::Real>& val, - const int start_ix, const int num_val) -{ - // call amrex::ParmParse::query, check if the user specified str. - std::vector<std::string> tmp_str_arr; - int is_specified = a_pp.queryarr(str, tmp_str_arr, start_ix, num_val); - if (is_specified) - { - // If so, create parser objects and apply them to the values provided by the user. - int const n = static_cast<int>(tmp_str_arr.size()); - val.resize(n); - for (int i=0 ; i < n ; i++) { - val[i] = static_cast<amrex::Real>(parseStringtoReal(tmp_str_arr[i])); - } - } - // return the same output as amrex::ParmParse::query - return is_specified; -} - -void -getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<amrex::Real>& val, - const int start_ix, const int num_val) -{ - // Create parser objects and apply them to the values provided by the user. - std::vector<std::string> tmp_str_arr; - a_pp.getarr(str, tmp_str_arr, start_ix, num_val); - - int const n = static_cast<int>(tmp_str_arr.size()); - val.resize(n); - for (int i=0 ; i < n ; i++) { - val[i] = static_cast<amrex::Real>(parseStringtoReal(tmp_str_arr[i])); - } -} - -int queryWithParser (const amrex::ParmParse& a_pp, char const * const str, int& val) { - amrex::Real rval; - const int result = queryWithParser(a_pp, str, rval); - if (result) { - val = safeCastToInt(std::round(rval), str); - } - return result; -} - -void getWithParser (const amrex::ParmParse& a_pp, char const * const str, int& val) { - amrex::Real rval; - getWithParser(a_pp, str, rval); - val = safeCastToInt(std::round(rval), str); -} - -int queryArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<int>& val, - const int start_ix, const int num_val) { - std::vector<amrex::Real> rval; - const int result = queryArrWithParser(a_pp, str, rval, start_ix, num_val); - if (result) { - val.resize(rval.size()); - for (unsigned long i = 0 ; i < val.size() ; i++) { - val[i] = safeCastToInt(std::round(rval[i]), str); - } - } - return result; -} - -void getArrWithParser (const amrex::ParmParse& a_pp, char const * const str, std::vector<int>& val, - const int start_ix, const int num_val) { - std::vector<amrex::Real> rval; - getArrWithParser(a_pp, str, rval, start_ix, num_val); - val.resize(rval.size()); - for (unsigned long i = 0 ; i < val.size() ; i++) { - val[i] = safeCastToInt(std::round(rval[i]), str); - } -} - void CheckDims () { // Ensure that geometry.dims is set properly. |