diff options
23 files changed, 349 insertions, 256 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 412299138..cfc9bae5e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -86,6 +86,13 @@ if(NOT WarpX_PRECISION IN_LIST WarpX_PRECISION_VALUES) message(FATAL_ERROR "WarpX_PRECISION (${WarpX_PRECISION}) must be one of ${WarpX_PRECISION_VALUES}") endif() +set(WarpX_PARTICLE_PRECISION_VALUES SINGLE DOUBLE) +set(WarpX_PARTICLE_PRECISION ${WarpX_PRECISION} CACHE STRING "Particle floating point precision (SINGLE/DOUBLE)") +set_property(CACHE WarpX_PARTICLE_PRECISION PROPERTY STRINGS ${WarpX_PARTICLE_PRECISION_VALUES}) +if(NOT WarpX_PARTICLE_PRECISION IN_LIST WarpX_PARTICLE_PRECISION_VALUES) + message(FATAL_ERROR "WarpX_PARTICLE_PRECISION (${WarpX_PARTICLE_PRECISION}) must be one of ${WarpX_PARTICLE_PRECISION_VALUES}") +endif() + set(WarpX_COMPUTE_VALUES NOACC OMP CUDA SYCL HIP) set(WarpX_COMPUTE OMP CACHE STRING "On-node, accelerated computing backend (NOACC/OMP/CUDA/SYCL/HIP)") set_property(CACHE WarpX_COMPUTE PROPERTY STRINGS ${WarpX_COMPUTE_VALUES}) diff --git a/Docs/source/install/cmake.rst b/Docs/source/install/cmake.rst index 22f997443..a50a05282 100644 --- a/Docs/source/install/cmake.rst +++ b/Docs/source/install/cmake.rst @@ -95,6 +95,7 @@ CMake Option Default & Values Descr ``WarpX_MPI_THREAD_MULTIPLE`` **ON**/OFF MPI thread-multiple support, i.e. for ``async_io`` ``WarpX_OPENPMD`` **ON**/OFF openPMD I/O (HDF5, ADIOS) ``WarpX_PRECISION`` SINGLE/**DOUBLE** Floating point precision (single/double) +``WarpX_PARTICLE_PRECISION`` SINGLE/**DOUBLE** Particle floating point precision (single/double), defaults to WarpX_PRECISION value if not set ``WarpX_PSATD`` ON/**OFF** Spectral solver ``WarpX_QED`` **ON**/OFF QED support (requires PICSAR) ``WarpX_QED_TABLE_GEN`` ON/**OFF** QED table generation support (requires PICSAR and Boost) @@ -266,6 +267,7 @@ Environment Variable Default & Values Descr ``WARPX_MPI`` ON/**OFF** Multi-node support (message-passing) ``WARPX_OPENPMD`` **ON**/OFF openPMD I/O (HDF5, ADIOS) ``WARPX_PRECISION`` SINGLE/**DOUBLE** Floating point precision (single/double) +``WARPX_PARTICLE_PRECISION`` SINGLE/**DOUBLE** Particle floating point precision (single/double), defaults to WarpX_PRECISION value if not set ``WARPX_PSATD`` ON/**OFF** Spectral solver ``WARPX_QED`` **ON**/OFF PICSAR QED (requires PICSAR) ``WARPX_QED_TABLE_GEN`` ON/**OFF** QED table generation (requires PICSAR and Boost) diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py index 629cecd2a..65a4d3cd5 100755 --- a/Python/pywarpx/_libwarpx.py +++ b/Python/pywarpx/_libwarpx.py @@ -133,19 +133,22 @@ class LibWarpX(): if _ParticleReal_size == 8: c_particlereal = ctypes.c_double - _numpy_particlereal_dtype = 'f8' + self._numpy_particlereal_dtype = 'f8' else: c_particlereal = ctypes.c_float - _numpy_particlereal_dtype = 'f4' + self._numpy_particlereal_dtype = 'f4' self.dim = self.libwarpx_so.warpx_SpaceDim() # our particle data type, depends on _ParticleReal_size - _p_struct = [(d, _numpy_particlereal_dtype) for d in 'xyz'[:self.dim]] + [('id', 'i4'), ('cpu', 'i4')] + _p_struct = ( + [(d, self._numpy_particlereal_dtype) for d in 'xyz'[:self.dim]] + + [('id', 'i4'), ('cpu', 'i4')] + ) self._p_dtype = np.dtype(_p_struct, align=True) _numpy_to_ctypes = {} - _numpy_to_ctypes[_numpy_particlereal_dtype] = c_particlereal + _numpy_to_ctypes[self._numpy_particlereal_dtype] = c_particlereal _numpy_to_ctypes['i4'] = ctypes.c_int class Particle(ctypes.Structure): @@ -594,32 +597,43 @@ class LibWarpX(): # --- Broadcast scalars into appropriate length arrays # --- If the parameter was not supplied, use the default value if lenx == 1: - x = np.full(maxlen, (x or 0.), float) + x = np.full(maxlen, (x or 0.), self._numpy_particlereal_dtype) if leny == 1: - y = np.full(maxlen, (y or 0.), float) + y = np.full(maxlen, (y or 0.), self._numpy_particlereal_dtype) if lenz == 1: - z = np.full(maxlen, (z or 0.), float) + z = np.full(maxlen, (z or 0.), self._numpy_particlereal_dtype) if lenux == 1: - ux = np.full(maxlen, (ux or 0.), float) + ux = np.full(maxlen, (ux or 0.), self._numpy_particlereal_dtype) if lenuy == 1: - uy = np.full(maxlen, (uy or 0.), float) + uy = np.full(maxlen, (uy or 0.), self._numpy_particlereal_dtype) if lenuz == 1: - uz = np.full(maxlen, (uz or 0.), float) + uz = np.full(maxlen, (uz or 0.), self._numpy_particlereal_dtype) if lenw == 1: - w = np.full(maxlen, (w or 0.), float) + w = np.full(maxlen, (w or 0.), self._numpy_particlereal_dtype) for key, val in kwargs.items(): if np.size(val) == 1: - kwargs[key] = np.full(maxlen, val, float) + kwargs[key] = np.full( + maxlen, val, self._numpy_particlereal_dtype + ) # --- The -3 is because the comps include the velocites nattr = self.get_nattr_species(species_name) - 3 - attr = np.zeros((maxlen, nattr)) + attr = np.zeros((maxlen, nattr), self._numpy_particlereal_dtype) attr[:,0] = w for key, vals in kwargs.items(): # --- The -3 is because components 1 to 3 are velocities attr[:,self.get_particle_comp_index(species_name, key)-3] = vals + # Iff x/y/z/ux/uy/uz are not numpy arrays of the correct dtype, new + # array copies are made with the correct dtype + x = x.astype(self._numpy_particlereal_dtype, copy=False) + y = y.astype(self._numpy_particlereal_dtype, copy=False) + z = z.astype(self._numpy_particlereal_dtype, copy=False) + ux = ux.astype(self._numpy_particlereal_dtype, copy=False) + uy = uy.astype(self._numpy_particlereal_dtype, copy=False) + uz = uz.astype(self._numpy_particlereal_dtype, copy=False) + self.libwarpx_so.warpx_addNParticles( ctypes.c_char_p(species_name.encode('utf-8')), x.size, x, y, z, ux, uy, uz, nattr, attr, unique_particles diff --git a/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json b/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json new file mode 100644 index 000000000..6e574d584 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/background_mcc_dp_psp.json @@ -0,0 +1,26 @@ +{ + "electrons": { + "particle_cpu": 187329.0, + "particle_id": 58191177610.0, + "particle_momentum_x": 1.014117400902589e-18, + "particle_momentum_y": 2.809684791375906e-19, + "particle_momentum_z": 2.802552971463722e-19, + "particle_position_x": 17126.5950836977, + "particle_position_y": 933.5069071217099, + "particle_weight": 61099470167.71875 + }, + "he_ions": { + "particle_cpu": 262776.0, + "particle_id": 206721076019.0, + "particle_momentum_x": 2.8874754507147993e-18, + "particle_momentum_y": 2.1946650983059055e-18, + "particle_momentum_z": 2.1983734218197867e-18, + "particle_position_x": 17604.642770509952, + "particle_position_y": 1099.9749558673466, + "particle_weight": 71967842052.5625 + }, + "lev=0": { + "rho_electrons": 0.035652566614896256, + "rho_he_ions": 0.04192071932755148 + } +}
\ No newline at end of file diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 3f0700fc9..0133c4c7e 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -3030,6 +3030,24 @@ compareParticles = 1 particleTypes = electrons he_ions analysisRoutine = Examples/analysis_default_regression.py +[background_mcc_dp_psp] +buildDir = . +inputFile = Examples/Physics_applications/capacitive_discharge/inputs_2d +runtime_params = warpx.abort_on_warning_threshold = high +dim = 2 +addToCompileString = +cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PRECISION=DOUBLE -DWarpX_PARTICLE_PRECISION=SINGLE -DWarpX_QED=OFF +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 1 +particleTypes = electrons he_ions +analysisRoutine = Examples/analysis_default_regression.py + [Python_background_mcc] buildDir = . inputFile = Examples/Physics_applications/capacitive_discharge/PICMI_inputs_2d.py 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. diff --git a/cmake/WarpXFunctions.cmake b/cmake/WarpXFunctions.cmake index c43d71d04..e57032b9e 100644 --- a/cmake/WarpXFunctions.cmake +++ b/cmake/WarpXFunctions.cmake @@ -208,6 +208,12 @@ function(set_warpx_binary_name) set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".SP") endif() + if(WarpX_PARTICLE_PRECISION STREQUAL "DOUBLE") + set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".PDP") + else() + set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".PSP") + endif() + if(WarpX_ASCENT) set_property(TARGET ${tgt} APPEND_STRING PROPERTY OUTPUT_NAME ".ASCENT") endif() @@ -374,6 +380,7 @@ function(warpx_print_summary) endif() message(" PSATD: ${WarpX_PSATD}") message(" PRECISION: ${WarpX_PRECISION}") + message(" PARTICLE PRECISION: ${WarpX_PARTICLE_PRECISION}") message(" OPENPMD: ${WarpX_OPENPMD}") message(" QED: ${WarpX_QED}") message(" QED table generation: ${WarpX_QED_TABLE_GEN}") diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 5aa791d64..80070d6f8 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -55,9 +55,13 @@ macro(find_amrex) if(WarpX_PRECISION STREQUAL "DOUBLE") set(AMReX_PRECISION "DOUBLE" CACHE INTERNAL "") - set(AMReX_PARTICLES_PRECISION "DOUBLE" CACHE INTERNAL "") else() set(AMReX_PRECISION "SINGLE" CACHE INTERNAL "") + endif() + + if(WarpX_PARTICLE_PRECISION STREQUAL "DOUBLE") + set(AMReX_PARTICLES_PRECISION "DOUBLE" CACHE INTERNAL "") + else() set(AMReX_PARTICLES_PRECISION "SINGLE" CACHE INTERNAL "") endif() @@ -221,7 +225,7 @@ macro(find_amrex) else() set(COMPONENT_SENSEI) endif() - set(COMPONENT_PRECISION ${WarpX_PRECISION} P${WarpX_PRECISION}) + set(COMPONENT_PRECISION ${WarpX_PRECISION} P${WarpX_PARTICLE_PRECISION}) find_package(AMReX 22.05 CONFIG REQUIRED COMPONENTS ${COMPONENT_ASCENT} ${COMPONENT_DIM} ${COMPONENT_EB} PARTICLES ${COMPONENT_PIC} ${COMPONENT_PRECISION} ${COMPONENT_SENSEI} TINYP LSOLVERS) message(STATUS "AMReX: Found version '${AMReX_VERSION}'") @@ -101,6 +101,7 @@ class CMakeBuild(build_ext): '-DWarpX_EB:BOOL=' + WARPX_EB, '-DWarpX_OPENPMD:BOOL=' + WARPX_OPENPMD, '-DWarpX_PRECISION=' + WARPX_PRECISION, + '-DWarpX_PARTICLE_PRECISION=' + WARPX_PARTICLE_PRECISION, '-DWarpX_PSATD:BOOL=' + WARPX_PSATD, '-DWarpX_QED:BOOL=' + WARPX_QED, '-DWarpX_QED_TABLE_GEN:BOOL=' + WARPX_QED_TABLE_GEN, @@ -201,6 +202,7 @@ WARPX_MPI = env.pop('WARPX_MPI', 'OFF') WARPX_EB = env.pop('WARPX_EB', 'OFF') WARPX_OPENPMD = env.pop('WARPX_OPENPMD', 'ON') WARPX_PRECISION = env.pop('WARPX_PRECISION', 'DOUBLE') +WARPX_PARTICLE_PRECISION = env.pop('WARPX_PARTICLE_PRECISION', WARPX_PRECISION) WARPX_PSATD = env.pop('WARPX_PSATD', 'OFF') WARPX_QED = env.pop('WARPX_QED', 'ON') WARPX_QED_TABLE_GEN = env.pop('WARPX_QED_TABLE_GEN', 'OFF') |