diff options
17 files changed, 313 insertions, 312 deletions
diff --git a/Source/Particles/Algorithms/KineticEnergy.H b/Source/Particles/Algorithms/KineticEnergy.H index af058178d..c14ff14fd 100644 --- a/Source/Particles/Algorithms/KineticEnergy.H +++ b/Source/Particles/Algorithms/KineticEnergy.H @@ -21,7 +21,7 @@ namespace Algorithms{ // This marks the gamma threshold to switch between the full relativistic expression // for particle kinetic energy and a Taylor expansion. static constexpr auto gamma_relativistic_threshold = - static_cast<amrex::Real>(1.005); + static_cast<amrex::ParticleReal>(1.005); /** * \brief Computes the kinetic energy of a particle. Below a threshold for the @@ -36,22 +36,22 @@ namespace Algorithms{ * @return the kinetic energy of the particle (in S.I. units) */ AMREX_GPU_HOST_DEVICE AMREX_INLINE - amrex::Real KineticEnergy( - const amrex::Real ux, const amrex::Real uy, const amrex::Real uz, - const amrex::Real mass) + amrex::ParticleReal KineticEnergy( + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, + const amrex::ParticleReal mass) { using namespace amrex; constexpr auto c2 = PhysConst::c * PhysConst::c; - constexpr auto inv_c2 = 1.0_rt/c2; + constexpr auto inv_c2 = 1.0_prt/c2; const auto u2 = (ux*ux + uy*uy + uz*uz)*inv_c2; - const auto gamma = std::sqrt(1.0_rt + u2); + const auto gamma = std::sqrt(1.0_prt + u2); const auto kk = (gamma > gamma_relativistic_threshold)? - (gamma-1.0_rt): - (u2*0.5_rt - u2*u2*(1.0_rt/8.0_rt) + u2*u2*u2*(1.0_rt/16.0_rt)- - u2*u2*u2*u2*(5.0_rt/128.0_rt) + (7.0_rt/256_rt)*u2*u2*u2*u2*u2); //Taylor expansion + (gamma-1.0_prt): + (u2*0.5_prt - u2*u2*(1.0_prt/8.0_prt) + u2*u2*u2*(1.0_prt/16.0_prt)- + u2*u2*u2*u2*(5.0_prt/128.0_prt) + (7.0_prt/256_prt)*u2*u2*u2*u2*u2); //Taylor expansion return kk*mass*c2; } @@ -66,8 +66,8 @@ namespace Algorithms{ * @return the kinetic energy of the photon (in S.I. units) */ AMREX_GPU_HOST_DEVICE AMREX_INLINE - amrex::Real KineticEnergyPhotons( - const amrex::Real ux, const amrex::Real uy, const amrex::Real uz) + amrex::ParticleReal KineticEnergyPhotons( + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz) { // Photons have zero mass, but ux, uy and uz are calculated assuming a mass equal to the // electron mass. Hence, photons need a special treatment to calculate the total energy. diff --git a/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.H b/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.H index e4ad27d04..d47234d45 100644 --- a/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.H +++ b/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.H @@ -27,7 +27,7 @@ public: virtual ~BackgroundMCCCollision () = default; - amrex::Real get_nu_max (amrex::Vector<MCCProcess> const& mcc_processes); + amrex::ParticleReal get_nu_max (amrex::Vector<MCCProcess> const& mcc_processes); /** Perform the collisions * @@ -73,14 +73,14 @@ private: bool init_flag = false; bool ionization_flag = false; - amrex::Real m_mass1; + amrex::ParticleReal m_mass1; - amrex::Real m_max_background_density = 0; - amrex::Real m_background_mass; - amrex::Real m_total_collision_prob; - amrex::Real m_total_collision_prob_ioniz = 0; - amrex::Real m_nu_max; - amrex::Real m_nu_max_ioniz; + amrex::ParticleReal m_max_background_density = 0; + amrex::ParticleReal m_background_mass; + amrex::ParticleReal m_total_collision_prob; + amrex::ParticleReal m_total_collision_prob_ioniz = 0; + amrex::ParticleReal m_nu_max; + amrex::ParticleReal m_nu_max_ioniz; amrex::Parser m_background_density_parser; amrex::Parser m_background_temperature_parser; diff --git a/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp b/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp index f9bfd9a93..1f45062b7 100644 --- a/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp +++ b/Source/Particles/Collision/BackgroundMCC/BackgroundMCCCollision.cpp @@ -28,7 +28,7 @@ BackgroundMCCCollision::BackgroundMCCCollision (std::string const collision_name amrex::ParmParse pp_collision_name(collision_name); - amrex::Real background_density = 0; + amrex::ParticleReal background_density = 0; if (queryWithParser(pp_collision_name, "background_density", background_density)) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (background_density > 0), "The background density must be greater than 0." @@ -41,7 +41,7 @@ BackgroundMCCCollision::BackgroundMCCCollision (std::string const collision_name m_background_density_parser = makeParser(background_density_str, {"x", "y", "z", "t"}); } - amrex::Real background_temperature; + amrex::ParticleReal background_temperature; if (queryWithParser(pp_collision_name, "background_temperature", background_temperature)) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE( (background_temperature >= 0), "The background temperature must be positive." @@ -89,7 +89,7 @@ BackgroundMCCCollision::BackgroundMCCCollision (std::string const collision_name std::string cross_section_file; pp_collision_name.query(kw_cross_section.c_str(), cross_section_file); - amrex::Real energy = 0.0; + amrex::ParticleReal energy = 0.0; // if the scattering process is excitation or ionization get the // energy associated with that process if (scattering_process.find("excitation") != std::string::npos || @@ -152,14 +152,14 @@ BackgroundMCCCollision::BackgroundMCCCollision (std::string const collision_name /** Calculate the maximum collision frequency using a fixed energy grid that * ranges from 1e-4 to 5000 eV in 0.2 eV increments */ -amrex::Real +amrex::ParticleReal BackgroundMCCCollision::get_nu_max(amrex::Vector<MCCProcess> const& mcc_processes) { using namespace amrex::literals; - amrex::Real nu, nu_max = 0.0; - amrex::Real E_start = 1e-4_rt; - amrex::Real E_end = 5000._rt; - amrex::Real E_step = 0.2_rt; + amrex::ParticleReal nu, nu_max = 0.0; + amrex::ParticleReal E_start = 1e-4_prt; + amrex::ParticleReal E_end = 5000._prt; + amrex::ParticleReal E_step = 0.2_prt; // set the energy limits and step size for calculating nu_max based // on the given cross-section inputs @@ -172,8 +172,8 @@ BackgroundMCCCollision::get_nu_max(amrex::Vector<MCCProcess> const& mcc_processe E_step = (energy_step < E_step) ? energy_step : E_step; } - for (amrex::Real E = E_start; E < E_end; E+=E_step) { - amrex::Real sigma_E = 0.0; + for (amrex::ParticleReal E = E_start; E < E_end; E+=E_step) { + amrex::ParticleReal sigma_E = 0.0; // loop through all collision pathways for (const auto &scattering_process : mcc_processes) { @@ -184,7 +184,7 @@ BackgroundMCCCollision::get_nu_max(amrex::Vector<MCCProcess> const& mcc_processe // calculate collision frequency nu = ( m_max_background_density - * std::sqrt(2.0_rt / m_mass1 * PhysConst::q_e) + * std::sqrt(2.0_prt / m_mass1 * PhysConst::q_e) * sigma_E * std::sqrt(E) ); if (nu > nu_max) { @@ -217,12 +217,12 @@ BackgroundMCCCollision::doCollisions (amrex::Real cur_time, amrex::Real dt, Mult // calculate total collision probability auto coll_n = m_nu_max * dt; - m_total_collision_prob = 1.0_rt - std::exp(-coll_n); + m_total_collision_prob = 1.0_prt - std::exp(-coll_n); // dt has to be small enough that a linear expansion of the collision // probability is sufficiently accurately, otherwise the MCC results // will be very heavily affected by small changes in the timestep - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n < 0.1_rt, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n < 0.1_prt, "dt is too large to ensure accurate MCC results" ); @@ -232,9 +232,9 @@ BackgroundMCCCollision::doCollisions (amrex::Real cur_time, amrex::Real dt, Mult // calculate total ionization probability auto coll_n_ioniz = m_nu_max_ioniz * dt; - m_total_collision_prob_ioniz = 1.0_rt - std::exp(-coll_n_ioniz); + m_total_collision_prob_ioniz = 1.0_prt - std::exp(-coll_n_ioniz); - WARPX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n_ioniz < 0.1_rt, + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n_ioniz < 0.1_prt, "dt is too large to ensure accurate MCC results" ); @@ -345,20 +345,20 @@ void BackgroundMCCCollision::doBackgroundCollisionsWithinTile amrex::ParticleReal x, y, z; GetPosition.AsStored(ip, x, y, z); - amrex::Real n_a = n_a_func(x, y, z, t); - amrex::Real T_a = T_a_func(x, y, z, t); + amrex::ParticleReal n_a = n_a_func(x, y, z, t); + amrex::ParticleReal T_a = T_a_func(x, y, z, t); - amrex::Real v_coll, v_coll2, sigma_E, nu_i = 0; + amrex::ParticleReal v_coll, v_coll2, sigma_E, nu_i = 0; double gamma, E_coll; amrex::ParticleReal ua_x, ua_y, ua_z, vx, vy, vz; amrex::ParticleReal uCOM_x, uCOM_y, uCOM_z; - amrex::Real col_select = amrex::Random(engine); + amrex::ParticleReal col_select = amrex::Random(engine); // get velocities of gas particles from a Maxwellian distribution auto const vel_std = sqrt(PhysConst::kb * T_a / M); - ua_x = vel_std * amrex::RandomNormal(0_rt, 1.0_rt, engine); - ua_y = vel_std * amrex::RandomNormal(0_rt, 1.0_rt, engine); - ua_z = vel_std * amrex::RandomNormal(0_rt, 1.0_rt, engine); + ua_x = vel_std * amrex::RandomNormal(0_prt, 1.0_prt, engine); + ua_y = vel_std * amrex::RandomNormal(0_prt, 1.0_prt, engine); + ua_z = vel_std * amrex::RandomNormal(0_prt, 1.0_prt, engine); // we assume the target particle is not relativistic (in // the lab frame) and therefore we can transform the projectile @@ -408,7 +408,7 @@ void BackgroundMCCCollision::doBackgroundCollisionsWithinTile // subtract any energy penalty of the collision from the // projectile energy - if (scattering_process.m_energy_penalty > 0.0_rt) { + if (scattering_process.m_energy_penalty > 0.0_prt) { ParticleUtils::getEnergy(v_coll2, m, E_coll); E_coll = (E_coll - scattering_process.m_energy_penalty) * PhysConst::q_e; auto scale_fac = sqrt(E_coll * (E_coll + 2.0_prt*mc2) / c2) / m / v_coll; @@ -464,7 +464,7 @@ void BackgroundMCCCollision::doBackgroundIonization m_nu_max_ioniz, m_background_density_func, t ); - amrex::Real sqrt_kb_m = std::sqrt(PhysConst::kb / m_background_mass); + amrex::ParticleReal sqrt_kb_m = std::sqrt(PhysConst::kb / m_background_mass); #ifdef AMREX_USE_OMP #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) diff --git a/Source/Particles/Collision/BackgroundMCC/ImpactIonization.H b/Source/Particles/Collision/BackgroundMCC/ImpactIonization.H index 8e6f1b049..0b48ae6b3 100644 --- a/Source/Particles/Collision/BackgroundMCC/ImpactIonization.H +++ b/Source/Particles/Collision/BackgroundMCC/ImpactIonization.H @@ -50,8 +50,8 @@ public: ImpactIonizationFilterFunc( MCCProcess const& mcc_process, double const mass, - amrex::Real const total_collision_prob, - amrex::Real const nu_max, + amrex::ParticleReal const total_collision_prob, + amrex::ParticleReal const nu_max, amrex::ParserExecutor<4> const& n_a_func, amrex::Real t ) : m_mcc_process(mcc_process.executor()), m_mass(mass), @@ -86,7 +86,7 @@ public: get_particle_position(p, x, y, z); // calculate neutral density at particle location - const Real n_a = m_n_a_func(x, y, z, m_t); + const ParticleReal n_a = m_n_a_func(x, y, z, m_t); // get the particle velocity const ParticleReal ux = ptd.m_rdata[PIdx::ux][i]; @@ -98,10 +98,10 @@ public: ParticleUtils::getEnergy(u_coll2, m_mass, E_coll); // get collision cross-section - const Real sigma_E = m_mcc_process.getCrossSection(E_coll); + const ParticleReal sigma_E = m_mcc_process.getCrossSection(E_coll); // calculate normalized collision frequency - const Real nu_i = n_a * sigma_E * sqrt(u_coll2) / m_nu_max; + const ParticleReal nu_i = n_a * sigma_E * sqrt(u_coll2) / m_nu_max; // check if this collision should be performed return (Random(engine) <= nu_i); @@ -110,8 +110,8 @@ public: private: MCCProcess::Executor m_mcc_process; double m_mass; - amrex::Real m_total_collision_prob = 0; - amrex::Real m_nu_max; + amrex::ParticleReal m_total_collision_prob = 0; + amrex::ParticleReal m_nu_max; amrex::ParserExecutor<4> m_n_a_func; amrex::Real m_t; }; @@ -146,7 +146,7 @@ public: * @param[in] t the current simulation time */ ImpactIonizationTransformFunc( - amrex::Real energy_cost, double mass1, amrex::Real sqrt_kb_m, + amrex::ParticleReal energy_cost, double mass1, amrex::ParticleReal sqrt_kb_m, amrex::ParserExecutor<4> const& T_a_func, amrex::Real t ) : m_energy_cost(energy_cost), m_mass1(mass1), m_sqrt_kb_m(sqrt_kb_m), m_T_a_func(T_a_func), m_t(t) { } @@ -182,7 +182,7 @@ public: // calculate standard deviation in neutral velocity distribution using // the local temperature - const Real ion_vel_std = m_sqrt_kb_m * sqrt(m_T_a_func(x, y, z, m_t)); + const ParticleReal ion_vel_std = m_sqrt_kb_m * sqrt(m_T_a_func(x, y, z, m_t)); // get references to the original particle's velocity auto& ux = src.m_rdata[PIdx::ux][i_src]; @@ -202,13 +202,13 @@ public: ParticleUtils::getEnergy(u_coll2, m_mass1, E_coll); // each electron gets half the energy (could change this later) - amrex::Real E_out = (E_coll - m_energy_cost) / 2.0_rt * PhysConst::q_e; + amrex::ParticleReal E_out = (E_coll - m_energy_cost) / 2.0_prt * PhysConst::q_e; // precalculate often used value constexpr auto c2 = PhysConst::c * PhysConst::c; const auto mc2 = m_mass1*c2; - amrex::Real up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m_mass1; + amrex::ParticleReal up = sqrt(E_out * (E_out + 2.0_prt*mc2) / c2) / m_mass1; // isotropically scatter electrons ParticleUtils::RandomizeVelocity(ux, uy, uz, up, engine); @@ -221,9 +221,9 @@ public: } private: - amrex::Real m_energy_cost; + amrex::ParticleReal m_energy_cost; double m_mass1; - amrex::Real m_sqrt_kb_m; + amrex::ParticleReal m_sqrt_kb_m; amrex::ParserExecutor<4> m_T_a_func; amrex::Real m_t; }; diff --git a/Source/Particles/Collision/BackgroundMCC/MCCProcess.H b/Source/Particles/Collision/BackgroundMCC/MCCProcess.H index 68ab61aba..0b28adf75 100644 --- a/Source/Particles/Collision/BackgroundMCC/MCCProcess.H +++ b/Source/Particles/Collision/BackgroundMCC/MCCProcess.H @@ -27,7 +27,7 @@ public: MCCProcess ( const std::string& scattering_process, const std::string& cross_section_file, - const amrex::Real energy + const amrex::ParticleReal energy ); template <typename InputVector> @@ -35,7 +35,7 @@ public: const std::string& scattering_process, const InputVector&& energies, const InputVector&& sigmas, - const amrex::Real energy + const amrex::ParticleReal energy ); MCCProcess (MCCProcess const&) = delete; @@ -55,14 +55,14 @@ public: static void readCrossSectionFile ( const std::string cross_section_file, - amrex::Vector<amrex::Real>& energies, - amrex::Gpu::HostVector<amrex::Real>& sigmas + amrex::Vector<amrex::ParticleReal>& energies, + amrex::Gpu::HostVector<amrex::ParticleReal>& sigmas ); static void sanityCheckEnergyGrid ( - const amrex::Vector<amrex::Real>& energies, - amrex::Real dE + const amrex::Vector<amrex::ParticleReal>& energies, + amrex::ParticleReal dE ); struct Executor { @@ -74,7 +74,7 @@ public: * */ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - amrex::Real getCrossSection (amrex::Real E_coll) const + amrex::ParticleReal getCrossSection (amrex::ParticleReal E_coll) const { if (E_coll < m_energy_lo) { return m_sigma_lo; @@ -84,7 +84,7 @@ public: using amrex::Math::floor; using amrex::Math::ceil; // calculate index of bounding energy pairs - amrex::Real temp = (E_coll - m_energy_lo) / m_dE; + amrex::ParticleReal temp = (E_coll - m_energy_lo) / m_dE; int idx_1 = static_cast<int>(floor(temp)); int idx_2 = static_cast<int>(ceil(temp)); @@ -94,9 +94,9 @@ public: } } - amrex::Real* m_sigmas_data = nullptr; - amrex::Real m_energy_lo, m_energy_hi, m_sigma_lo, m_sigma_hi, m_dE; - amrex::Real m_energy_penalty; + amrex::ParticleReal* m_sigmas_data = nullptr; + amrex::ParticleReal m_energy_lo, m_energy_hi, m_sigma_lo, m_sigma_hi, m_dE; + amrex::ParticleReal m_energy_penalty; MCCProcessType m_type; }; @@ -108,15 +108,15 @@ public: #endif } - amrex::Real getCrossSection (amrex::Real E_coll) const + amrex::ParticleReal getCrossSection (amrex::ParticleReal E_coll) const { return m_exe_h.getCrossSection(E_coll); } - amrex::Real getEnergyPenalty () const { return m_exe_h.m_energy_penalty; } - amrex::Real getMinEnergyInput () const { return m_exe_h.m_energy_lo; } - amrex::Real getMaxEnergyInput () const { return m_exe_h.m_energy_hi; } - amrex::Real getEnergyInputStep () const { return m_exe_h.m_dE; } + amrex::ParticleReal getEnergyPenalty () const { return m_exe_h.m_energy_penalty; } + amrex::ParticleReal getMinEnergyInput () const { return m_exe_h.m_energy_lo; } + amrex::ParticleReal getMaxEnergyInput () const { return m_exe_h.m_energy_hi; } + amrex::ParticleReal getEnergyInputStep () const { return m_exe_h.m_dE; } MCCProcessType type () const { return m_exe_h.m_type; } @@ -125,15 +125,15 @@ private: static MCCProcessType parseProcessType(const std::string& process); - void init (const std::string& scattering_process, const amrex::Real energy); + void init (const std::string& scattering_process, const amrex::ParticleReal energy); - amrex::Vector<amrex::Real> m_energies; + amrex::Vector<amrex::ParticleReal> m_energies; #ifdef AMREX_USE_GPU - amrex::Gpu::DeviceVector<amrex::Real> m_sigmas_d; + amrex::Gpu::DeviceVector<amrex::ParticleReal> m_sigmas_d; Executor m_exe_d; #endif - amrex::Gpu::HostVector<amrex::Real> m_sigmas_h; + amrex::Gpu::HostVector<amrex::ParticleReal> m_sigmas_h; Executor m_exe_h; int m_grid_size; diff --git a/Source/Particles/Collision/BackgroundMCC/MCCProcess.cpp b/Source/Particles/Collision/BackgroundMCC/MCCProcess.cpp index b9e6f12e6..0f6dae3bc 100644 --- a/Source/Particles/Collision/BackgroundMCC/MCCProcess.cpp +++ b/Source/Particles/Collision/BackgroundMCC/MCCProcess.cpp @@ -12,7 +12,7 @@ MCCProcess::MCCProcess ( const std::string& scattering_process, const std::string& cross_section_file, - const amrex::Real energy ) + const amrex::ParticleReal energy ) { // read the cross-section data file into memory readCrossSectionFile(cross_section_file, m_energies, m_sigmas_h); @@ -25,7 +25,7 @@ MCCProcess::MCCProcess ( const std::string& scattering_process, const InputVector&& energies, const InputVector&& sigmas, - const amrex::Real energy ) + const amrex::ParticleReal energy ) { m_energies.insert(m_energies.begin(), std::begin(energies), std::end(energies)); m_sigmas_h.insert(m_sigmas_h.begin(), std::begin(sigmas), std::end(sigmas)); @@ -34,7 +34,7 @@ MCCProcess::MCCProcess ( } void -MCCProcess::init (const std::string& scattering_process, const amrex::Real energy) +MCCProcess::init (const std::string& scattering_process, const amrex::ParticleReal energy) { using namespace amrex::literals; m_exe_h.m_sigmas_data = m_sigmas_h.data(); @@ -45,7 +45,7 @@ MCCProcess::init (const std::string& scattering_process, const amrex::Real energ m_exe_h.m_energy_hi = m_energies[m_grid_size-1]; m_exe_h.m_sigma_lo = m_sigmas_h[0]; m_exe_h.m_sigma_hi = m_sigmas_h[m_grid_size-1]; - m_exe_h.m_dE = (m_exe_h.m_energy_hi - m_exe_h.m_energy_lo)/(m_grid_size - 1._rt); + m_exe_h.m_dE = (m_exe_h.m_energy_hi - m_exe_h.m_energy_lo)/(m_grid_size - 1._prt); m_exe_h.m_energy_penalty = energy; m_exe_h.m_type = parseProcessType(scattering_process); @@ -93,13 +93,13 @@ MCCProcess::parseProcessType(const std::string& scattering_process) void MCCProcess::readCrossSectionFile ( const std::string cross_section_file, - amrex::Vector<amrex::Real>& energies, - amrex::Gpu::HostVector<amrex::Real>& sigmas ) + amrex::Vector<amrex::ParticleReal>& energies, + amrex::Gpu::HostVector<amrex::ParticleReal>& sigmas ) { std::ifstream infile(cross_section_file); if(!infile.is_open()) amrex::Abort("Failed to open cross-section data file"); - amrex::Real energy, sigma; + amrex::ParticleReal energy, sigma; while (infile >> energy >> sigma) { energies.push_back(energy); sigmas.push_back(sigma); @@ -110,8 +110,8 @@ MCCProcess::readCrossSectionFile ( void MCCProcess::sanityCheckEnergyGrid ( - const amrex::Vector<amrex::Real>& energies, - amrex::Real dE + const amrex::Vector<amrex::ParticleReal>& energies, + amrex::ParticleReal dE ) { // confirm that the input data for the cross-section was provided with diff --git a/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H index 56538102f..ee4b3b016 100644 --- a/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H +++ b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.H @@ -47,7 +47,7 @@ public: * */ void doBackgroundStoppingOnElectronsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t, - amrex::Real species_mass, amrex::Real species_charge); + amrex::ParticleReal species_mass, amrex::ParticleReal species_charge); /** Perform the stopping calculation within a tile for stopping on ions * @@ -59,12 +59,12 @@ public: * */ void doBackgroundStoppingOnIonsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t, - amrex::Real species_mass, amrex::Real species_charge); + amrex::ParticleReal species_mass, amrex::ParticleReal species_charge); private: - amrex::Real m_background_mass; - amrex::Real m_background_charge_state; + amrex::ParticleReal m_background_mass; + amrex::ParticleReal m_background_charge_state; BackgroundStoppingType m_background_type; amrex::Parser m_background_density_parser; diff --git a/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp index a970c490f..939f098ee 100644 --- a/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp +++ b/Source/Particles/Collision/BackgroundStopping/BackgroundStopping.cpp @@ -33,7 +33,7 @@ BackgroundStopping::BackgroundStopping (std::string const collision_name) AMREX_ALWAYS_ASSERT_WITH_MESSAGE(false, "background_type must be either electrons or ions"); } - amrex::Real background_density; + amrex::ParticleReal background_density; std::string background_density_str; if (queryWithParser(pp_collision_name, "background_density", background_density)) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(background_density > 0, @@ -46,7 +46,7 @@ BackgroundStopping::BackgroundStopping (std::string const collision_name) "For background stopping, the background density must be specified."); } - amrex::Real background_temperature; + amrex::ParticleReal background_temperature; std::string background_temperature_str; if (queryWithParser(pp_collision_name, "background_temperature", background_temperature)) { AMREX_ALWAYS_ASSERT_WITH_MESSAGE(background_temperature > 0, @@ -81,8 +81,8 @@ BackgroundStopping::doCollisions (amrex::Real cur_time, amrex::Real dt, MultiPar using namespace amrex::literals; auto& species = mypc->GetParticleContainerFromName(m_species_names[0]); - amrex::Real species_mass = species.getMass(); - amrex::Real species_charge = species.getCharge(); + amrex::ParticleReal species_mass = species.getMass(); + amrex::ParticleReal species_charge = species.getCharge(); BackgroundStoppingType background_type = m_background_type; @@ -121,7 +121,7 @@ BackgroundStopping::doCollisions (amrex::Real cur_time, amrex::Real dt, MultiPar } void BackgroundStopping::doBackgroundStoppingOnElectronsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t, - amrex::Real species_mass, amrex::Real species_charge) + amrex::ParticleReal species_mass, amrex::ParticleReal species_charge) { using namespace amrex::literals; @@ -132,7 +132,7 @@ void BackgroundStopping::doBackgroundStoppingOnElectronsWithinTile (WarpXParIter const long np = pti.numParticles(); // get background particle mass - amrex::Real mass_e = m_background_mass; + amrex::ParticleReal mass_e = m_background_mass; // setup parsers for the background density and temperature auto n_e_func = m_background_density_func; @@ -153,33 +153,33 @@ void BackgroundStopping::doBackgroundStoppingOnElectronsWithinTile (WarpXParIter amrex::ParticleReal x, y, z; GetPosition.AsStored(ip, x, y, z); - amrex::Real const n_e = n_e_func(x, y, z, t); - amrex::Real const T_e = T_e_func(x, y, z, t)*PhysConst::kb; + amrex::ParticleReal const n_e = n_e_func(x, y, z, t); + amrex::ParticleReal const T_e = T_e_func(x, y, z, t)*PhysConst::kb; // This implements the equation 14.12 from Introduction to Plasma Physics, // Goldston and Rutherford, the slowing down of beam ions due to collisions with electrons. // The equation is written as dV/dt = -alpha*V, and integrated to // give V(t+dt) = V(t)*exp(-alpha*dt) - amrex::Real constexpr pi = MathConst::pi; - amrex::Real constexpr ep0 = PhysConst::ep0; - amrex::Real constexpr q_e = PhysConst::q_e; - amrex::Real constexpr q_e2 = q_e*q_e; - amrex::Real constexpr ep02 = ep0*ep0; + amrex::ParticleReal constexpr pi = MathConst::pi; + amrex::ParticleReal constexpr ep0 = PhysConst::ep0; + amrex::ParticleReal constexpr q_e = PhysConst::q_e; + amrex::ParticleReal constexpr q_e2 = q_e*q_e; + amrex::ParticleReal constexpr ep02 = ep0*ep0; - amrex::Real const Zb = abs(species_charge/q_e); + amrex::ParticleReal const Zb = abs(species_charge/q_e); - amrex::Real const vth = sqrt(3._rt*T_e/mass_e); - amrex::Real const wp = sqrt(n_e*q_e2/(ep0*mass_e)); - amrex::Real const lambdadb = vth/wp; - amrex::Real const lambdadb3 = lambdadb*lambdadb*lambdadb; - amrex::Real const loglambda = log((12._rt*pi/Zb)*(n_e*lambdadb3)); + amrex::ParticleReal const vth = sqrt(3._prt*T_e/mass_e); + amrex::ParticleReal const wp = sqrt(n_e*q_e2/(ep0*mass_e)); + amrex::ParticleReal const lambdadb = vth/wp; + amrex::ParticleReal const lambdadb3 = lambdadb*lambdadb*lambdadb; + amrex::ParticleReal const loglambda = log((12._prt*pi/Zb)*(n_e*lambdadb3)); - amrex::Real const pi32 = pi*sqrt(pi); - amrex::Real const q2 = species_charge*species_charge; - amrex::Real const T32 = T_e*sqrt(T_e); + amrex::ParticleReal const pi32 = pi*sqrt(pi); + amrex::ParticleReal const q2 = species_charge*species_charge; + amrex::ParticleReal const T32 = T_e*sqrt(T_e); - amrex::Real const alpha = sqrt(2._rt)*n_e*q2*q_e2*sqrt(mass_e)*loglambda/(12._rt*pi32*ep02*species_mass*T32); + amrex::ParticleReal const alpha = sqrt(2._prt)*n_e*q2*q_e2*sqrt(mass_e)*loglambda/(12._prt*pi32*ep02*species_mass*T32); ux[ip] *= exp(-alpha*dt); uy[ip] *= exp(-alpha*dt); @@ -190,7 +190,7 @@ void BackgroundStopping::doBackgroundStoppingOnElectronsWithinTile (WarpXParIter } void BackgroundStopping::doBackgroundStoppingOnIonsWithinTile (WarpXParIter& pti, amrex::Real dt, amrex::Real t, - amrex::Real species_mass, amrex::Real species_charge) + amrex::ParticleReal species_mass, amrex::ParticleReal species_charge) { using namespace amrex::literals; @@ -201,8 +201,8 @@ void BackgroundStopping::doBackgroundStoppingOnIonsWithinTile (WarpXParIter& pti const long np = pti.numParticles(); // get background particle mass - amrex::Real mass_i = m_background_mass; - amrex::Real charge_state_i = m_background_charge_state; + amrex::ParticleReal mass_i = m_background_mass; + amrex::ParticleReal charge_state_i = m_background_charge_state; // setup parsers for the background density and temperature auto n_i_func = m_background_density_func; @@ -223,37 +223,37 @@ void BackgroundStopping::doBackgroundStoppingOnIonsWithinTile (WarpXParIter& pti amrex::ParticleReal x, y, z; GetPosition.AsStored(ip, x, y, z); - amrex::Real const n_i = n_i_func(x, y, z, t); - amrex::Real const T_i = T_i_func(x, y, z, t)*PhysConst::kb; + amrex::ParticleReal const n_i = n_i_func(x, y, z, t); + amrex::ParticleReal const T_i = T_i_func(x, y, z, t)*PhysConst::kb; // This implements the equation 14.20 from Introduction to Plasma Physics, // Goldston and Rutherford, the slowing down of beam ions due to collisions with electrons. // The equation is written with energy, W, as dW/dt = -alpha/W**0.5, and integrated to // give W(t+dt) = (W(t)**1.5 - 3./2.*alpha*dt)**(2/3) - amrex::Real constexpr pi = MathConst::pi; - amrex::Real constexpr q_e = PhysConst::q_e; - amrex::Real constexpr q_e2 = q_e*q_e; - amrex::Real constexpr ep0 = PhysConst::ep0; - amrex::Real constexpr ep02 = ep0*ep0; + amrex::ParticleReal constexpr pi = MathConst::pi; + amrex::ParticleReal constexpr q_e = PhysConst::q_e; + amrex::ParticleReal constexpr q_e2 = q_e*q_e; + amrex::ParticleReal constexpr ep0 = PhysConst::ep0; + amrex::ParticleReal constexpr ep02 = ep0*ep0; - amrex::Real const qi2 = charge_state_i*charge_state_i*q_e2; - amrex::Real const qb2 = species_charge*species_charge; - amrex::Real const Zb = abs(species_charge/q_e); + amrex::ParticleReal const qi2 = charge_state_i*charge_state_i*q_e2; + amrex::ParticleReal const qb2 = species_charge*species_charge; + amrex::ParticleReal const Zb = abs(species_charge/q_e); - amrex::Real const vth = sqrt(3._rt*T_i/mass_i); - amrex::Real const wp = sqrt(n_i*q_e2/(ep0*mass_i)); - amrex::Real const lambdadb = vth/wp; - amrex::Real const lambdadb3 = lambdadb*lambdadb*lambdadb; - amrex::Real const loglambda = log((12._rt*pi/Zb)*(n_i*lambdadb3)); + amrex::ParticleReal const vth = sqrt(3._prt*T_i/mass_i); + amrex::ParticleReal const wp = sqrt(n_i*q_e2/(ep0*mass_i)); + amrex::ParticleReal const lambdadb = vth/wp; + amrex::ParticleReal const lambdadb3 = lambdadb*lambdadb*lambdadb; + amrex::ParticleReal const loglambda = log((12._prt*pi/Zb)*(n_i*lambdadb3)); - amrex::Real const alpha = sqrt(2._rt)*n_i*qi2*qb2*sqrt(species_mass)*loglambda/(8._rt*pi*ep02*mass_i); + amrex::ParticleReal const alpha = sqrt(2._prt)*n_i*qi2*qb2*sqrt(species_mass)*loglambda/(8._prt*pi*ep02*mass_i); - amrex::Real const W0 = 0.5_rt*species_mass*(ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip]); - amrex::Real const f1 = pow(W0, 1.5_rt) - 1.5_rt*alpha*dt; + amrex::ParticleReal const W0 = 0.5_prt*species_mass*(ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip]); + amrex::ParticleReal const f1 = pow(W0, 1.5_prt) - 1.5_prt*alpha*dt; // If f1 goes negative, the particle has fully stopped, so set W1 to 0. - amrex::Real const W1 = pow((f1 > 0._rt ? f1 : 0._rt), 2._rt/3._rt); - amrex::Real const vscale = (W0 > 0._rt ? std::sqrt(W1/W0) : 0._rt); + amrex::ParticleReal const W1 = pow((f1 > 0._prt ? f1 : 0._prt), 2._prt/3._prt); + amrex::ParticleReal const vscale = (W0 > 0._prt ? std::sqrt(W1/W0) : 0._prt); ux[ip] *= vscale; uy[ip] *= vscale; diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index b1c1ebd93..0c5ec182e 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -339,7 +339,7 @@ public: indices_1, cell_start_1, cell_half_1, engine ); #if defined WARPX_DIM_RZ int ri = (i_cell - i_cell%nz) / nz; - auto dV = MathConst::pi*(2.0_rt*ri+1.0_rt)*dr*dr*dz; + auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz; #endif // Call the function in order to perform collisions // If there are product species, mask, p_pair_indices_1/2, and @@ -500,7 +500,7 @@ public: ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine); #if defined WARPX_DIM_RZ int ri = (i_cell - i_cell%nz) / nz; - auto dV = MathConst::pi*(2.0_rt*ri+1.0_rt)*dr*dr*dz; + auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz; #endif // Call the function in order to perform collisions // If there are product species, p_mask, p_pair_indices_1/2, and diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H b/Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H index 1aafba313..c38eb6819 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/ElasticCollisionPerez.H @@ -18,7 +18,8 @@ /** Prepare information for and call UpdateMomentumPerezElastic(). * * @tparam T_index type of index arguments - * @tparam T_R type of floating point arguments + * @tparam T_PR type of particle related floating point arguments + * @tparam T_R type of other floating point arguments * @tparam SoaData_type type of the "struct of array" for the two involved species * @param[in] I1s,I2s is the start index for I1,I2 (inclusive). * @param[in] I1e,I2e is the stop index for I1,I2 (exclusive). @@ -37,7 +38,7 @@ * @param[in] engine the random number generator state & factory */ -template <typename T_index, typename T_R, typename SoaData_type> +template <typename T_index, typename T_PR, typename T_R, typename SoaData_type> AMREX_GPU_HOST_DEVICE AMREX_INLINE void ElasticCollisionPerez ( T_index const I1s, T_index const I1e, @@ -45,42 +46,42 @@ void ElasticCollisionPerez ( T_index const* AMREX_RESTRICT I1, T_index const* AMREX_RESTRICT I2, SoaData_type soa_1, SoaData_type soa_2, - T_R const q1, T_R const q2, - T_R const m1, T_R const m2, - T_R const T1, T_R const T2, - T_R const dt, T_R const L, T_R const dV, + T_PR const q1, T_PR const q2, + T_PR const m1, T_PR const m2, + T_PR const T1, T_PR const T2, + T_R const dt, T_PR const L, T_R const dV, amrex::RandomEngine const& engine) { int NI1 = I1e - I1s; int NI2 = I2e - I2s; - T_R * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; - T_R * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; - T_R * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; - T_R * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; + T_PR * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; + T_PR * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; + T_PR * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; + T_PR * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; - T_R * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; - T_R * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux]; - T_R * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy]; - T_R * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; + T_PR * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; + T_PR * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux]; + T_PR * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy]; + T_PR * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; // get local T1t and T2t - T_R T1t; T_R T2t; - if ( T1 <= T_R(0.0) && L <= T_R(0.0) ) + T_PR T1t; T_PR T2t; + if ( T1 <= T_PR(0.0) && L <= T_PR(0.0) ) { T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1); } else { T1t = T1; } - if ( T2 <= T_R(0.0) && L <= T_R(0.0) ) + if ( T2 <= T_PR(0.0) && L <= T_PR(0.0) ) { T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2); } else { T2t = T2; } // local density - T_R n1 = T_R(0.0); - T_R n2 = T_R(0.0); - T_R n12 = T_R(0.0); + T_PR n1 = T_PR(0.0); + T_PR n2 = T_PR(0.0); + T_PR n12 = T_PR(0.0); for (int i1=I1s; i1<static_cast<int>(I1e); ++i1) { n1 += w1[ I1[i1] ]; } for (int i2=I2s; i2<static_cast<int>(I2e); ++i2) { n2 += w2[ I2[i2] ]; } n1 = n1 / dV; n2 = n2 / dV; @@ -96,21 +97,21 @@ void ElasticCollisionPerez ( } // compute Debye length lmdD - T_R lmdD; - if ( T1t < T_R(0.0) || T2t < T_R(0.0) ) { - lmdD = T_R(0.0); + T_PR lmdD; + if ( T1t < T_PR(0.0) || T2t < T_PR(0.0) ) { + lmdD = T_PR(0.0); } else { - lmdD = T_R(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) + - n2*q2*q2/(T2t*PhysConst::ep0) ); + lmdD = T_PR(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) + + n2*q2*q2/(T2t*PhysConst::ep0) ); } - T_R rmin = std::pow( T_R(4.0) * MathConst::pi / T_R(3.0) * - amrex::max(n1,n2), T_R(-1.0/3.0) ); + T_PR rmin = std::pow( T_PR(4.0) * MathConst::pi / T_PR(3.0) * + amrex::max(n1,n2), T_PR(-1.0/3.0) ); lmdD = amrex::max(lmdD, rmin); #if (defined WARPX_DIM_RZ) - T_R * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; - T_R * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta]; + T_PR * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; + T_PR * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta]; #endif // call UpdateMomentumPerezElastic() @@ -128,8 +129,8 @@ void ElasticCollisionPerez ( * momentum of one of the macroparticles in agreement with this cylindrical symmetry. * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation; * there is a corresponding assert statement at initialization.) */ - T_R const theta = theta2[I2[i2]]-theta1[I1[i1]]; - T_R const u1xbuf = u1x[I1[i1]]; + T_PR const theta = theta2[I2[i2]]-theta1[I1[i1]]; + T_PR const u1xbuf = u1x[I1[i1]]; u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta); u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta); #endif @@ -143,7 +144,7 @@ void ElasticCollisionPerez ( engine); #if (defined WARPX_DIM_RZ) - T_R const u1xbuf_new = u1x[I1[i1]]; + T_PR const u1xbuf_new = u1x[I1[i1]]; u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta); u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta); #endif diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index 3387f606a..0acdc9af8 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -48,7 +48,7 @@ public: using namespace amrex::literals; amrex::ParmParse pp_collision_name(collision_name); // default Coulomb log, if < 0, will be computed automatically - m_CoulombLog = -1.0_rt; + m_CoulombLog = -1.0_prt; queryWithParser(pp_collision_name, "CoulombLog", m_CoulombLog); } @@ -76,7 +76,7 @@ public: GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, amrex::ParticleReal const q1, amrex::ParticleReal const q2, amrex::ParticleReal const m1, amrex::ParticleReal const m2, - amrex::ParticleReal const dt, amrex::ParticleReal const dV, + amrex::Real const dt, amrex::Real 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*/, diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H index fb7dc4da1..c24736ec0 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/UpdateMomentumPerezElastic.H @@ -25,60 +25,60 @@ * compile with USE_ASSERTION=TRUE. */ -template <typename T_R> +template <typename T_PR, typename T_R> AMREX_GPU_HOST_DEVICE AMREX_INLINE void UpdateMomentumPerezElastic ( - T_R& u1x, T_R& u1y, T_R& u1z, T_R& u2x, T_R& u2y, T_R& u2z, - T_R const n1, T_R const n2, T_R const n12, - T_R const q1, T_R const m1, T_R const w1, - T_R const q2, T_R const m2, T_R const w2, - T_R const dt, T_R const L, T_R const lmdD, + T_PR& u1x, T_PR& u1y, T_PR& u1z, T_PR& u2x, T_PR& u2y, T_PR& u2z, + T_PR const n1, T_PR const n2, T_PR const n12, + T_PR const q1, T_PR const m1, T_PR const w1, + T_PR const q2, T_PR const m2, T_PR const w2, + T_R const dt, T_PR const L, T_PR const lmdD, amrex::RandomEngine const& engine) { - T_R const diffx = amrex::Math::abs(u1x-u2x); - T_R const diffy = amrex::Math::abs(u1y-u2y); - T_R const diffz = amrex::Math::abs(u1z-u2z); - T_R const diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz); - T_R const summm = std::sqrt(u1x*u1x+u1y*u1y+u1z*u1z) + std::sqrt(u2x*u2x+u2y*u2y+u2z*u2z); + T_PR const diffx = amrex::Math::abs(u1x-u2x); + T_PR const diffy = amrex::Math::abs(u1y-u2y); + T_PR const diffz = amrex::Math::abs(u1z-u2z); + T_PR const diffm = std::sqrt(diffx*diffx+diffy*diffy+diffz*diffz); + T_PR const summm = std::sqrt(u1x*u1x+u1y*u1y+u1z*u1z) + std::sqrt(u2x*u2x+u2y*u2y+u2z*u2z); // If g = u1 - u2 = 0, do not collide. // Or if the relative difference is less than 1.0e-10. - if ( diffm < std::numeric_limits<T_R>::min() || diffm/summm < 1.0e-10 ) { return; } + if ( diffm < std::numeric_limits<T_PR>::min() || diffm/summm < 1.0e-10 ) { return; } - T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c ); + T_PR constexpr inv_c2 = T_PR(1.0) / ( PhysConst::c * PhysConst::c ); // Compute Lorentz factor gamma - T_R const g1 = std::sqrt( T_R(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 ); - T_R const g2 = std::sqrt( T_R(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 ); + T_PR const g1 = std::sqrt( T_PR(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 ); + T_PR const g2 = std::sqrt( T_PR(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 ); // Compute momenta - T_R const p1x = u1x * m1; - T_R const p1y = u1y * m1; - T_R const p1z = u1z * m1; - T_R const p2x = u2x * m2; - T_R const p2y = u2y * m2; - T_R const p2z = u2z * m2; + T_PR const p1x = u1x * m1; + T_PR const p1y = u1y * m1; + T_PR const p1z = u1z * m1; + T_PR const p2x = u2x * m2; + T_PR const p2y = u2y * m2; + T_PR const p2z = u2z * m2; // Compute center-of-mass (COM) velocity and gamma - T_R const mass_g = m1 * g1 + m2 * g2; - T_R const vcx = (p1x+p2x) / mass_g; - T_R const vcy = (p1y+p2y) / mass_g; - T_R const vcz = (p1z+p2z) / mass_g; - T_R const vcms = vcx*vcx + vcy*vcy + vcz*vcz; - T_R const gc = T_R(1.0) / std::sqrt( T_R(1.0) - vcms*inv_c2 ); + T_PR const mass_g = m1 * g1 + m2 * g2; + T_PR const vcx = (p1x+p2x) / mass_g; + T_PR const vcy = (p1y+p2y) / mass_g; + T_PR const vcz = (p1z+p2z) / mass_g; + T_PR const vcms = vcx*vcx + vcy*vcy + vcz*vcz; + T_PR const gc = T_PR(1.0) / std::sqrt( T_PR(1.0) - vcms*inv_c2 ); // Compute vc dot v1 and v2 - T_R const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1; - T_R const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2; + T_PR const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1; + T_PR const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2; // Compute p1 star - T_R p1sx; - T_R p1sy; - T_R p1sz; - if ( vcms > std::numeric_limits<T_R>::min() ) + T_PR p1sx; + T_PR p1sy; + T_PR p1sz; + if ( vcms > std::numeric_limits<T_PR>::min() ) { - T_R const lorentz_tansform_factor = - ( (gc-T_R(1.0))/vcms*vcDv1 - gc )*m1*g1; + T_PR const lorentz_tansform_factor = + ( (gc-T_PR(1.0))/vcms*vcDv1 - gc )*m1*g1; p1sx = p1x + vcx*lorentz_tansform_factor; p1sy = p1y + vcy*lorentz_tansform_factor; p1sz = p1z + vcz*lorentz_tansform_factor; @@ -89,45 +89,45 @@ void UpdateMomentumPerezElastic ( p1sy = p1y; p1sz = p1z; } - T_R const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz ); + T_PR const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz ); // Compute gamma star - T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1; - T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2; + T_PR const g1s = ( T_PR(1.0) - vcDv1*inv_c2 )*gc*g1; + T_PR const g2s = ( T_PR(1.0) - vcDv2*inv_c2 )*gc*g2; // Compute the Coulomb log lnLmd - T_R lnLmd; - if ( L > T_R(0.0) ) { lnLmd = L; } + T_PR lnLmd; + if ( L > T_PR(0.0) ) { lnLmd = L; } else { // Compute b0 - T_R const b0 = amrex::Math::abs(q1*q2) * inv_c2 / - (T_R(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g * - ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_R(1.0) ); + T_PR const b0 = amrex::Math::abs(q1*q2) * inv_c2 / + (T_PR(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g * + ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_PR(1.0) ); // Compute the minimal impact parameter - constexpr T_R hbar_pi = static_cast<T_R>(PhysConst::hbar*MathConst::pi); - const T_R bmin = amrex::max(hbar_pi/p1sm, b0); + constexpr T_PR hbar_pi = static_cast<T_PR>(PhysConst::hbar*MathConst::pi); + const T_PR bmin = amrex::max(hbar_pi/p1sm, b0); // Compute the Coulomb log lnLmd - lnLmd = amrex::max( T_R(2.0), - T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) ); + lnLmd = amrex::max( T_PR(2.0), + T_PR(0.5)*std::log(T_PR(1.0)+lmdD*lmdD/(bmin*bmin)) ); } // Compute s - const auto tts = m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_R(1.0); + const auto tts = m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_PR(1.0); const auto tts2 = tts*tts; - T_R s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / - ( T_R(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * + T_PR s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / + ( T_PR(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * tts2; // Compute s' const auto cbrt_n1 = std::cbrt(n1); const auto cbrt_n2 = std::cbrt(n2); - const auto coeff = static_cast<T_R>( + const auto coeff = static_cast<T_PR>( std::pow(4.0*MathConst::pi/3.0,1.0/3.0)); - T_R const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc); - T_R const sp = coeff * n1*n2/n12 * dt * vrel * (m1+m2) / + T_PR const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc); + T_PR const sp = coeff * n1*n2/n12 * dt * vrel * (m1+m2) / amrex::max( m1*cbrt_n1*cbrt_n1, m2*cbrt_n2*cbrt_n2); @@ -135,54 +135,54 @@ void UpdateMomentumPerezElastic ( s = amrex::min(s,sp); // Get random numbers - T_R r = amrex::Random(engine); + T_PR r = amrex::Random(engine); // Compute scattering angle - T_R cosXs; - T_R sinXs; - if ( s <= T_R(0.1) ) + T_PR cosXs; + T_PR sinXs; + if ( s <= T_PR(0.1) ) { while ( true ) { - cosXs = T_R(1.0) + s * std::log(r); + cosXs = T_PR(1.0) + s * std::log(r); // Avoid the bug when r is too small such that cosXs < -1 - if ( cosXs >= T_R(-1.0) ) { break; } + if ( cosXs >= T_PR(-1.0) ) { break; } r = amrex::Random(engine); } } - else if ( s > T_R(0.1) && s <= T_R(3.0) ) + else if ( s > T_PR(0.1) && s <= T_PR(3.0) ) { - T_R const Ainv = static_cast<T_R>( + T_PR const Ainv = static_cast<T_PR>( 0.0056958 + 0.9560202*s - 0.508139*s*s + 0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s); - cosXs = Ainv * std::log( std::exp(T_R(-1.0)/Ainv) + - T_R(2.0) * r * std::sinh(T_R(1.0)/Ainv) ); + cosXs = Ainv * std::log( std::exp(T_PR(-1.0)/Ainv) + + T_PR(2.0) * r * std::sinh(T_PR(1.0)/Ainv) ); } - else if ( s > T_R(3.0) && s <= T_R(6.0) ) + else if ( s > T_PR(3.0) && s <= T_PR(6.0) ) { - T_R const A = T_R(3.0) * std::exp(-s); - cosXs = T_R(1.0)/A * std::log( std::exp(-A) + - T_R(2.0) * r * std::sinh(A) ); + T_PR const A = T_PR(3.0) * std::exp(-s); + cosXs = T_PR(1.0)/A * std::log( std::exp(-A) + + T_PR(2.0) * r * std::sinh(A) ); } else { - cosXs = T_R(2.0) * r - T_R(1.0); + cosXs = T_PR(2.0) * r - T_PR(1.0); } - sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs); + sinXs = std::sqrt(T_PR(1.0) - cosXs*cosXs); // Get random azimuthal angle - T_R const phis = amrex::Random(engine) * T_R(2.0) * MathConst::pi; - T_R const cosphis = std::cos(phis); - T_R const sinphis = std::sin(phis); + T_PR const phis = amrex::Random(engine) * T_PR(2.0) * MathConst::pi; + T_PR const cosphis = std::cos(phis); + T_PR const sinphis = std::sin(phis); // Compute post-collision momenta pfs in COM - T_R p1fsx; - T_R p1fsy; - T_R p1fsz; + T_PR p1fsx; + T_PR p1fsy; + T_PR p1fsz; // p1sp is the p1s perpendicular - T_R p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy ); + T_PR p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy ); // Make sure p1sp is not almost zero - if ( p1sp > std::numeric_limits<T_R>::min() ) + if ( p1sp > std::numeric_limits<T_PR>::min() ) { p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis + ( p1sy*p1sm/p1sp ) * sinXs*sinphis + @@ -191,7 +191,7 @@ void UpdateMomentumPerezElastic ( (-p1sx*p1sm/p1sp ) * sinXs*sinphis + ( p1sy ) * cosXs; p1fsz = (-p1sp ) * sinXs*cosphis + - ( T_R(0.0) ) * sinXs*sinphis + + ( T_PR(0.0) ) * sinXs*sinphis + ( p1sz ) * cosXs; // Note a negative sign is different from // Eq. (12) in Perez's paper, @@ -210,25 +210,25 @@ void UpdateMomentumPerezElastic ( (-p1sy*p1sm/p1sp ) * sinXs*sinphis + ( p1sz ) * cosXs; p1fsx = (-p1sp ) * sinXs*cosphis + - ( T_R(0.0) ) * sinXs*sinphis + + ( T_PR(0.0) ) * sinXs*sinphis + ( p1sx ) * cosXs; } - T_R const p2fsx = -p1fsx; - T_R const p2fsy = -p1fsy; - T_R const p2fsz = -p1fsz; + T_PR const p2fsx = -p1fsx; + T_PR const p2fsy = -p1fsy; + T_PR const p2fsz = -p1fsz; // Transform from COM to lab frame - T_R p1fx; T_R p2fx; - T_R p1fy; T_R p2fy; - T_R p1fz; T_R p2fz; - if ( vcms > std::numeric_limits<T_R>::min() ) + T_PR p1fx; T_PR p2fx; + T_PR p1fy; T_PR p2fy; + T_PR p1fz; T_PR p2fz; + if ( vcms > std::numeric_limits<T_PR>::min() ) { - T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz; - T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz; - T_R const factor = (gc-T_R(1.0))/vcms; - T_R const factor1 = factor*vcDp1fs + m1*g1s*gc; - T_R const factor2 = factor*vcDp2fs + m2*g2s*gc; + T_PR const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz; + T_PR const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz; + T_PR const factor = (gc-T_PR(1.0))/vcms; + T_PR const factor1 = factor*vcDp1fs + m1*g1s*gc; + T_PR const factor2 = factor*vcDp2fs + m2*g2s*gc; p1fx = p1fsx + vcx * factor1; p1fy = p1fsy + vcy * factor1; p1fz = p1fsz + vcz * factor1; diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index 2da560ea3..a807a4262 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -93,14 +93,14 @@ public: } // default fusion multiplier - m_fusion_multiplier = 1.0_rt; + m_fusion_multiplier = 1.0_prt; queryWithParser(pp_collision_name, "fusion_multiplier", m_fusion_multiplier); // default fusion probability threshold - m_probability_threshold = 0.02_rt; + m_probability_threshold = 0.02_prt; queryWithParser(pp_collision_name, "fusion_probability_threshold", m_probability_threshold); // default fusion probability target_value - m_probability_target_value = 0.002_rt; + m_probability_target_value = 0.002_prt; queryWithParser(pp_collision_name, "fusion_probability_target_value", m_probability_target_value); } @@ -223,15 +223,15 @@ public: private: // Factor used to increase the number of fusion reaction by decreasing the weight of the // produced particles - amrex::Real m_fusion_multiplier; + amrex::ParticleReal m_fusion_multiplier; // If the fusion multiplier is too high and results in a fusion probability that approaches // 1, there is a risk of underestimating the total fusion yield. In these cases, we reduce // the fusion multiplier used in a given collision. m_probability_threshold is the fusion // probability threshold above which we reduce the fusion multiplier. // m_probability_target_value is the target probability used to determine by how much // the fusion multiplier should be reduced. - amrex::Real m_probability_threshold; - amrex::Real m_probability_target_value; + amrex::ParticleReal m_probability_threshold; + amrex::ParticleReal m_probability_target_value; NuclearFusionType m_fusion_type; bool m_isSameSpecies; }; diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H index 4df0c583e..0bab4bef4 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H @@ -56,7 +56,7 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part const amrex::ParticleReal& u2y, const amrex::ParticleReal& u2z, const amrex::ParticleReal& m1, const amrex::ParticleReal& m2, amrex::ParticleReal w1, amrex::ParticleReal w2, - const amrex::Real& dt, const amrex::Real& dV, const int& pair_index, + const amrex::Real& dt, const amrex::ParticleReal& dV, const int& pair_index, index_type* AMREX_RESTRICT p_mask, amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, const amrex::ParticleReal& fusion_multiplier, diff --git a/Source/Particles/Gather/GetExternalFields.H b/Source/Particles/Gather/GetExternalFields.H index 66b983853..c7123c9f5 100644 --- a/Source/Particles/Gather/GetExternalFields.H +++ b/Source/Particles/Gather/GetExternalFields.H @@ -28,8 +28,8 @@ struct GetExternalEBField ExternalFieldInitType m_Etype; ExternalFieldInitType m_Btype; - amrex::Real m_gamma_boost; - amrex::Real m_uz_boost; + amrex::ParticleReal m_gamma_boost; + amrex::ParticleReal m_uz_boost; amrex::GpuArray<amrex::ParticleReal, 3> m_Efield_value; amrex::GpuArray<amrex::ParticleReal, 3> m_Bfield_value; @@ -44,11 +44,11 @@ struct GetExternalEBField GetParticlePosition m_get_position; amrex::Real m_time; - amrex::Real m_repeated_plasma_lens_period; - const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_starts = nullptr; - const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_lengths = nullptr; - const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_strengths_E = nullptr; - const amrex::Real* AMREX_RESTRICT m_repeated_plasma_lens_strengths_B = nullptr; + amrex::ParticleReal m_repeated_plasma_lens_period; + const amrex::ParticleReal* AMREX_RESTRICT m_repeated_plasma_lens_starts = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_repeated_plasma_lens_lengths = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_repeated_plasma_lens_strengths_E = nullptr; + const amrex::ParticleReal* AMREX_RESTRICT m_repeated_plasma_lens_strengths_B = nullptr; int m_n_lenses; amrex::Real m_dt; const amrex::ParticleReal* AMREX_RESTRICT m_ux = nullptr; @@ -68,14 +68,14 @@ struct GetExternalEBField if (m_Etype == None && m_Btype == None) return; - amrex::ParticleReal Ex = 0._rt; - amrex::ParticleReal Ey = 0._rt; - amrex::ParticleReal Ez = 0._rt; - amrex::ParticleReal Bx = 0._rt; - amrex::ParticleReal By = 0._rt; - amrex::ParticleReal Bz = 0._rt; + amrex::ParticleReal Ex = 0._prt; + amrex::ParticleReal Ey = 0._prt; + amrex::ParticleReal Ez = 0._prt; + amrex::ParticleReal Bx = 0._prt; + amrex::ParticleReal By = 0._prt; + amrex::ParticleReal Bz = 0._prt; - constexpr amrex::Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c); + constexpr amrex::ParticleReal inv_c2 = 1._prt/(PhysConst::c*PhysConst::c); if (m_Etype == Constant) { @@ -88,13 +88,13 @@ struct GetExternalEBField amrex::ParticleReal x, y, z; m_get_position(i, x, y, z); amrex::Real lab_time = m_time; - if (m_gamma_boost > 1._rt) { + if (m_gamma_boost > 1._prt) { 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((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); + Ex = m_Exfield_partparser((amrex::ParticleReal) x, (amrex::ParticleReal) y, (amrex::ParticleReal) z, lab_time); + Ey = m_Eyfield_partparser((amrex::ParticleReal) x, (amrex::ParticleReal) y, (amrex::ParticleReal) z, lab_time); + Ez = m_Ezfield_partparser((amrex::ParticleReal) x, (amrex::ParticleReal) y, (amrex::ParticleReal) z, lab_time); } if (m_Btype == Constant) @@ -108,7 +108,7 @@ struct GetExternalEBField amrex::ParticleReal x, y, z; m_get_position(i, x, y, z); amrex::Real lab_time = m_time; - if (m_gamma_boost > 1._rt) { + if (m_gamma_boost > 1._prt) { lab_time = m_gamma_boost*m_time + m_uz_boost*z*inv_c2; z = m_gamma_boost*z + m_uz_boost*m_time; } @@ -127,13 +127,13 @@ struct GetExternalEBField const amrex::ParticleReal uyp = m_uy[i]; const amrex::ParticleReal uzp = m_uz[i]; - const amrex::ParticleReal gamma = std::sqrt(1._rt + (uxp*uxp + uyp*uyp + uzp*uzp)*inv_c2); + const amrex::ParticleReal gamma = std::sqrt(1._prt + (uxp*uxp + uyp*uyp + uzp*uzp)*inv_c2); const amrex::ParticleReal vzp = uzp/gamma; amrex::ParticleReal zl = z; amrex::ParticleReal zr = z + vzp*m_dt; - if (m_gamma_boost > 1._rt) { + if (m_gamma_boost > 1._prt) { zl = m_gamma_boost*zl + m_uz_boost*m_time; zr = m_gamma_boost*zr + m_uz_boost*(m_time + m_dt); } @@ -141,18 +141,18 @@ struct GetExternalEBField // This assumes that zl > 0. int i_lens = static_cast<int>(std::floor(zl/m_repeated_plasma_lens_period)); i_lens = i_lens % m_n_lenses; - amrex::Real const lens_start = m_repeated_plasma_lens_starts[i_lens] + i_lens*m_repeated_plasma_lens_period; - amrex::Real const lens_end = lens_start + m_repeated_plasma_lens_lengths[i_lens]; + amrex::ParticleReal const lens_start = m_repeated_plasma_lens_starts[i_lens] + i_lens*m_repeated_plasma_lens_period; + amrex::ParticleReal const lens_end = lens_start + m_repeated_plasma_lens_lengths[i_lens]; // Calculate the residence correction // frac will be 1 if the step is completely inside the lens, between 0 and 1 // when entering or leaving the lens, and otherwise 0. // This assumes that vzp > 0. - amrex::Real fl = 0.; + amrex::ParticleReal fl = 0.; if (zl >= lens_start && zl < lens_end) fl = 1.; - amrex::Real fr = 0.; + amrex::ParticleReal fr = 0.; if (zr >= lens_start && zr < lens_end) fr = 1.; - amrex::Real frac = fl; + amrex::ParticleReal frac = fl; if (fl > fr) frac = (lens_end - zl)/(zr - zl); if (fr > fl) frac = (zr - lens_start)/(zr - zl); @@ -165,12 +165,12 @@ struct GetExternalEBField } - if (m_gamma_boost > 1._rt) { + if (m_gamma_boost > 1._prt) { // Transform the fields to the boosted frame - const amrex::Real Ex_boost = m_gamma_boost*Ex - m_uz_boost*By; - const amrex::Real Ey_boost = m_gamma_boost*Ey + m_uz_boost*Bx; - const amrex::Real Bx_boost = m_gamma_boost*Bx + m_uz_boost*Ey*inv_c2; - const amrex::Real By_boost = m_gamma_boost*By - m_uz_boost*Ex*inv_c2; + const amrex::ParticleReal Ex_boost = m_gamma_boost*Ex - m_uz_boost*By; + const amrex::ParticleReal Ey_boost = m_gamma_boost*Ey + m_uz_boost*Bx; + const amrex::ParticleReal Bx_boost = m_gamma_boost*Bx + m_uz_boost*Ey*inv_c2; + const amrex::ParticleReal By_boost = m_gamma_boost*By - m_uz_boost*Ex*inv_c2; Ex = Ex_boost; Ey = Ey_boost; Bx = Bx_boost; diff --git a/Source/Particles/Gather/GetExternalFields.cpp b/Source/Particles/Gather/GetExternalFields.cpp index 8c57ab9cc..6f5c4cea2 100644 --- a/Source/Particles/Gather/GetExternalFields.cpp +++ b/Source/Particles/Gather/GetExternalFields.cpp @@ -17,7 +17,7 @@ GetExternalEBField::GetExternalEBField (const WarpXParIter& a_pti, int a_offset) auto& mypc = warpx.GetPartContainer(); m_gamma_boost = WarpX::gamma_boost; - m_uz_boost = std::sqrt(WarpX::gamma_boost*WarpX::gamma_boost - 1._rt)*PhysConst::c; + m_uz_boost = std::sqrt(WarpX::gamma_boost*WarpX::gamma_boost - 1._prt)*PhysConst::c; m_Etype = Unknown; m_Btype = Unknown; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 80a2f362d..88981c110 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -298,8 +298,8 @@ public: std::string m_B_ext_particle_s = "none"; std::string m_E_ext_particle_s = "none"; // External fields added to particle fields. - amrex::Vector<amrex::Real> m_B_external_particle; - amrex::Vector<amrex::Real> m_E_external_particle; + amrex::Vector<amrex::ParticleReal> m_B_external_particle; + amrex::Vector<amrex::ParticleReal> m_E_external_particle; // Parser for B_external on the particle std::unique_ptr<amrex::Parser> m_Bx_particle_parser; std::unique_ptr<amrex::Parser> m_By_particle_parser; @@ -309,15 +309,15 @@ public: std::unique_ptr<amrex::Parser> m_Ey_particle_parser; std::unique_ptr<amrex::Parser> m_Ez_particle_parser; - amrex::Real m_repeated_plasma_lens_period; - amrex::Vector<amrex::Real> h_repeated_plasma_lens_starts; - amrex::Vector<amrex::Real> h_repeated_plasma_lens_lengths; - amrex::Vector<amrex::Real> h_repeated_plasma_lens_strengths_E; - amrex::Vector<amrex::Real> h_repeated_plasma_lens_strengths_B; - amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_starts; - amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_lengths; - amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_strengths_E; - amrex::Gpu::DeviceVector<amrex::Real> d_repeated_plasma_lens_strengths_B; + amrex::ParticleReal m_repeated_plasma_lens_period; + amrex::Vector<amrex::ParticleReal> h_repeated_plasma_lens_starts; + amrex::Vector<amrex::ParticleReal> h_repeated_plasma_lens_lengths; + amrex::Vector<amrex::ParticleReal> h_repeated_plasma_lens_strengths_E; + amrex::Vector<amrex::ParticleReal> h_repeated_plasma_lens_strengths_B; + amrex::Gpu::DeviceVector<amrex::ParticleReal> d_repeated_plasma_lens_starts; + amrex::Gpu::DeviceVector<amrex::ParticleReal> d_repeated_plasma_lens_lengths; + amrex::Gpu::DeviceVector<amrex::ParticleReal> d_repeated_plasma_lens_strengths_E; + amrex::Gpu::DeviceVector<amrex::ParticleReal> d_repeated_plasma_lens_strengths_B; #ifdef WARPX_QED /** |