diff options
Diffstat (limited to 'Source/Particles/Collision')
| -rw-r--r-- | Source/Particles/Collision/BackgroundMCCCollision.H | 77 | ||||
| -rw-r--r-- | Source/Particles/Collision/BackgroundMCCCollision.cpp | 382 | ||||
| -rw-r--r-- | Source/Particles/Collision/CMakeLists.txt | 2 | ||||
| -rw-r--r-- | Source/Particles/Collision/CollisionHandler.cpp | 7 | ||||
| -rw-r--r-- | Source/Particles/Collision/MCCProcess.H | 116 | ||||
| -rw-r--r-- | Source/Particles/Collision/MCCProcess.cpp | 129 | ||||
| -rw-r--r-- | Source/Particles/Collision/MCCScattering.H | 307 | ||||
| -rw-r--r-- | Source/Particles/Collision/Make.package | 2 |
8 files changed, 1022 insertions, 0 deletions
diff --git a/Source/Particles/Collision/BackgroundMCCCollision.H b/Source/Particles/Collision/BackgroundMCCCollision.H new file mode 100644 index 000000000..c739667a9 --- /dev/null +++ b/Source/Particles/Collision/BackgroundMCCCollision.H @@ -0,0 +1,77 @@ +/* Copyright 2021 Modern Electron + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PARTICLES_COLLISION_BACKGROUNDMCCCOLLISION_H_ +#define WARPX_PARTICLES_COLLISION_BACKGROUNDMCCCOLLISION_H_ + +#include "Particles/MultiParticleContainer.H" +#include "CollisionBase.H" +#include "MCCProcess.H" + +#include <AMReX_REAL.H> +#include <AMReX_GpuContainers.H> + +#include <string> + +class BackgroundMCCCollision + : public CollisionBase +{ +public: + template <typename T> + using VectorType = amrex::Gpu::ManagedVector<T>; + + BackgroundMCCCollision (std::string collision_name); + + amrex::Real get_nu_max ( VectorType<MCCProcess*> const& mcc_processes ); + + /** Perform the collisions + * + * @param cur_time Current time + * @param mypc Container of species involved + * + */ + void doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc) override; + + /** Perform particle conserving MCC collisions within a tile + * + * @param pti particle iterator + * + */ + void doBackgroundCollisionsWithinTile ( WarpXParIter& pti ); + + /** Perform MCC ionization interactions + * + * @param pti particle iterator + * @param species1/2 reference to species container used to inject + * new particles + * + */ + void doBackgroundIonization ( + int lev, + WarpXParticleContainer& species1, + WarpXParticleContainer& species2 + ); + +private: + + VectorType< MCCProcess* > m_scattering_processes; + VectorType< MCCProcess* > m_ionization_processes; + + bool init_flag = false; + bool ionization_flag = false; + + amrex::Real m_mass1; + + amrex::Real m_background_temperature; + amrex::Real m_background_density; + 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; +}; + +#endif // WARPX_PARTICLES_COLLISION_BACKGROUNDMCCCOLLISION_H_ diff --git a/Source/Particles/Collision/BackgroundMCCCollision.cpp b/Source/Particles/Collision/BackgroundMCCCollision.cpp new file mode 100644 index 000000000..c723a2800 --- /dev/null +++ b/Source/Particles/Collision/BackgroundMCCCollision.cpp @@ -0,0 +1,382 @@ +/* Copyright 2021 Modern Electron + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "BackgroundMCCCollision.H" +#include "MCCScattering.H" +#include "Particles/ParticleCreation/FilterCopyTransform.H" +#include "Particles/ParticleCreation/SmartCopy.H" +#include "Utils/ParticleUtils.H" +#include "Utils/WarpXUtil.H" +#include "Utils/WarpXProfilerWrapper.H" +#include "WarpX.H" + +#include <AMReX_ParmParse.H> +#include <AMReX_REAL.H> +#include <AMReX_Vector.H> + +#include <string> + +BackgroundMCCCollision::BackgroundMCCCollision (std::string const collision_name) + : CollisionBase(collision_name) +{ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_species_names.size() == 1, + "Background MCC must have exactly one species."); + + amrex::ParmParse pp(collision_name); + + pp.query("background_density", m_background_density); + pp.query("background_temperature", m_background_temperature); + + // if the neutral mass is specified use it, but if ionization is + // included the mass of the secondary species of that interaction + // will be used. If no neutral mass is specified and ionization is not + // included the mass of the colliding species will be used + m_background_mass = -1; + pp.query("background_mass", m_background_mass); + + // query for a list of collision processes + // these could be elastic, excitation, charge_exchange, back, etc. + amrex::Vector<std::string> scattering_process_names; + pp.queryarr("scattering_processes", scattering_process_names); + + // create a vector of MCCProcess objects from each scattering + // process name + for (auto scattering_process : scattering_process_names) { + std::string kw_cross_section = scattering_process + "_cross_section"; + std::string cross_section_file; + pp.query(kw_cross_section.c_str(), cross_section_file); + + amrex::Real 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 || + scattering_process.find("ionization") != std::string::npos) { + std::string kw_energy = scattering_process + "_energy"; + pp.get(kw_energy.c_str(), energy); + } + + auto process = new MCCProcess(scattering_process, cross_section_file, energy); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(process->m_type != MCCProcessType::INVALID, + "Cannot add an unknown MCC process type"); + + // if the scattering process is ionization get the secondary species + // only one ionization process is supported, the vector + // m_ionization_processes is only used to make it simple to calculate + // the maximum collision frequency with the same function used for + // particle conserving processes + if (process->m_type == MCCProcessType::IONIZATION) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!ionization_flag, + "Background MCC only supports a single ionization process"); + ionization_flag = true; + + std::string secondary_species; + pp.get("ionization_species", secondary_species); + m_species_names.push_back(secondary_species); + + m_ionization_processes.push_back(process); + } else { + m_scattering_processes.push_back(process); + } + + amrex::Gpu::synchronize(); + } +} + +/** 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 +BackgroundMCCCollision::get_nu_max(amrex::Gpu::ManagedVector<MCCProcess*> const& mcc_processes) +{ + using namespace amrex::literals; + amrex::Real nu, nu_max = 0.0; + + for (double E = 1e-4; E < 5000; E+=0.2) { + amrex::Real sigma_E = 0.0; + + // loop through all collision pathways + for (const auto &scattering_process : mcc_processes) { + // get collision cross-section + sigma_E += scattering_process->getCrossSection(E); + } + + // calculate collision frequency + nu = ( + m_background_density * std::sqrt(2.0_rt / m_mass1 * PhysConst::q_e) + * sigma_E * std::sqrt(E) + ); + if (nu > nu_max) { + nu_max = nu; + } + } + return nu_max; +} + +void +BackgroundMCCCollision::doCollisions (amrex::Real cur_time, MultiParticleContainer* mypc) +{ + WARPX_PROFILE("BackgroundMCCCollision::doCollisions()"); + using namespace amrex::literals; + + const amrex::Real dt = WarpX::GetInstance().getdt(0); + if ( int(std::floor(cur_time/dt)) % m_ndt != 0 ) return; + + auto& species1 = mypc->GetParticleContainerFromName(m_species_names[0]); + // this is a very ugly hack to have species2 be a reference and be + // defined in the scope of doCollisions + auto& species2 = ( + (m_species_names.size() == 2) ? + mypc->GetParticleContainerFromName(m_species_names[1]) : + mypc->GetParticleContainerFromName(m_species_names[0]) + ); + + if (!init_flag) { + m_mass1 = species1.getMass(); + + // calculate maximum collision frequency without ionization + m_nu_max = get_nu_max(m_scattering_processes); + + // calculate total collision probability + auto coll_n = m_nu_max * dt; + m_total_collision_prob = 1.0_rt - 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 + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n < 0.1_rt, + "dt is too large to ensure accurate MCC results" + ); + + if (ionization_flag) { + // calculate maximum collision frequency for ionization + m_nu_max_ioniz = get_nu_max(m_ionization_processes); + + // 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); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(coll_n_ioniz < 0.1_rt, + "dt is too large to ensure accurate MCC results" + ); + + // if an ionization process is included the secondary species mass + // is taken as the background mass + m_background_mass = species2.getMass(); + } + // if no neutral species mass was specified and ionization is not + // included assume that the collisions will be with neutrals of the + // same mass as the colliding species (like with ion-neutral collisions) + else if (m_background_mass == -1) { + m_background_mass = species1.getMass(); + } + + amrex::Print() << + "Setting up collisions for " << m_species_names[0] << " with total " + "collision probability: " << + m_total_collision_prob + m_total_collision_prob_ioniz << "\n"; + + init_flag = true; + } + + // Loop over refinement levels + auto const flvl = species1.finestLevel(); + for (int lev = 0; lev <= flvl; ++lev) { + + // firstly loop over particles box by box and do all particle conserving + // scattering +#ifdef _OPENMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (WarpXParIter pti(species1, lev); pti.isValid(); ++pti) { + doBackgroundCollisionsWithinTile(pti); + } + + // secondly perform ionization through the SmartCopyFactory if needed + if (ionization_flag) { + doBackgroundIonization(lev, species1, species2); + } + } +} + + +/** Perform all particle conserving MCC collisions within a tile + * + * @param pti particle iterator + * + */ +void BackgroundMCCCollision::doBackgroundCollisionsWithinTile +( WarpXParIter& pti ) +{ + using namespace amrex::literals; + + // So that CUDA code gets its intrinsic, not the host-only C++ library version + using std::sqrt; + + // get particle count + const long np = pti.numParticles(); + + // get collider properties + amrex::Real mass1 = m_mass1; + + // get neutral properties + amrex::Real n_a = m_background_density; + amrex::Real T_a = m_background_temperature; + amrex::Real mass_a = m_background_mass; + amrex::Real vel_std = sqrt(PhysConst::kb * T_a / mass_a); + + // get collision parameters + auto scattering_processes = m_scattering_processes.data(); + auto process_count = m_scattering_processes.size(); + + amrex::Real total_collision_prob = m_total_collision_prob; + amrex::Real nu_max = m_nu_max; + + // get Struct-Of-Array particle data, also called attribs + auto& attribs = pti.GetAttribs(); + amrex::ParticleReal* const AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); + amrex::ParticleReal* const AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); + + amrex::ParallelForRNG(np, + [=] AMREX_GPU_HOST_DEVICE (long ip, amrex::RandomEngine const& engine) + { + // determine if this particle should collide + if (amrex::Random(engine) > total_collision_prob) return; + + amrex::Real v_coll, v_coll2, E_coll, sigma_E, nu_i = 0; + amrex::Real col_select = amrex::Random(engine); + amrex::ParticleReal ua_x, ua_y, ua_z; + amrex::ParticleReal uCOM_x, uCOM_y, uCOM_z; + + // get velocities of gas particles from a Maxwellian distribution + 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); + + // calculate the center of momentum velocity + uCOM_x = (mass1 * ux[ip] + mass_a * ua_x) / (mass1 + mass_a); + uCOM_y = (mass1 * uy[ip] + mass_a * ua_y) / (mass1 + mass_a); + uCOM_z = (mass1 * uz[ip] + mass_a * ua_z) / (mass1 + mass_a); + + // calculate relative velocity of collision and collision energy if + // the colliding particle is an ion. For electron collisions we + // cannot use the relative velocity since that allows the + // possibility where the electron kinetic energy in the lab frame + // is insufficient to cause excitation but not in the COM frame - + // for energy to balance this situation requires the neutral to + // lose energy during the collision which we don't currently + // account for. + if (mass_a / mass1 > 1e3) { + v_coll2 = ux[ip]*ux[ip] + uy[ip]*uy[ip] + uz[ip]*uz[ip]; + E_coll = 0.5_rt * mass1 * v_coll2 / PhysConst::q_e; + } + else { + v_coll2 = ( + (ux[ip] - ua_x)*(ux[ip] - ua_x) + + (uy[ip] - ua_y)*(uy[ip] - ua_y) + + (uz[ip] - ua_z)*(uz[ip] - ua_z) + ); + E_coll = ( + 0.5_rt * mass1 * mass_a / (mass1 + mass_a) * v_coll2 + / PhysConst::q_e + ); + } + v_coll = sqrt(v_coll2); + + // loop through all collision pathways + for (size_t i = 0; i < process_count; i++) { + auto const& scattering_process = **(scattering_processes + i); + + // get collision cross-section + sigma_E = scattering_process.getCrossSection(E_coll); + + // calculate normalized collision frequency + nu_i += n_a * sigma_E * v_coll / nu_max; + + // check if this collision should be performed and call + // the appropriate scattering function + if (col_select > nu_i) continue; + + if (scattering_process.m_type == MCCProcessType::ELASTIC) { + ElasticScattering( + ux[ip], uy[ip], uz[ip], uCOM_x, uCOM_y, uCOM_z, engine + ); + } + else if (scattering_process.m_type == MCCProcessType::BACK) { + BackScattering( + ux[ip], uy[ip], uz[ip], uCOM_x, uCOM_y, uCOM_z + ); + } + else if (scattering_process.m_type == MCCProcessType::CHARGE_EXCHANGE) { + ChargeExchange(ux[ip], uy[ip], uz[ip], ua_x, ua_y, ua_z); + } + else if (scattering_process.m_type == MCCProcessType::EXCITATION) { + // get the new velocity magnitude + amrex::Real vp = sqrt( + 2.0_rt / mass1 * PhysConst::q_e + * (E_coll - scattering_process.m_energy_penalty) + ); + RandomizeVelocity(ux[ip], uy[ip], uz[ip], vp, engine); + } + break; + } + } + ); +} + + +/** Perform MCC ionization interactions + * + * @param pti particle iterator + * @param species1/2 reference to species container used to inject new + particles from ionization events + * + */ +void BackgroundMCCCollision::doBackgroundIonization +( int lev, WarpXParticleContainer& species1, + WarpXParticleContainer& species2) +{ + WARPX_PROFILE("BackgroundMCCCollision::doBackgroundIonization()"); + + SmartCopyFactory copy_factory_elec(species1, species1); + SmartCopyFactory copy_factory_ion(species1, species2); + const auto CopyElec = copy_factory_elec.getSmartCopy(); + const auto CopyIon = copy_factory_ion.getSmartCopy(); + + const auto Filter = ImpactIonizationFilterFunc( + *m_ionization_processes[0], + m_mass1, m_total_collision_prob_ioniz, + m_nu_max_ioniz / m_background_density + ); + + amrex::Real vel_std = std::sqrt( + PhysConst::kb * m_background_temperature / m_background_mass + ); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + for (WarpXParIter pti(species1, lev); pti.isValid(); ++pti) { + auto& elec_tile = species1.ParticlesAt(lev, pti); + auto& ion_tile = species2.ParticlesAt(lev, pti); + + const auto np_elec = elec_tile.numParticles(); + const auto np_ion = ion_tile.numParticles(); + + auto Transform = ImpactIonizationTransformFunc( + m_ionization_processes[0]->m_energy_penalty, m_mass1, vel_std + ); + + const auto num_added = filterCopyTransformParticles<1>( + elec_tile, ion_tile, elec_tile, np_elec, np_ion, + Filter, CopyElec, CopyIon, Transform + ); + + setNewParticleIDs(elec_tile, np_elec, num_added); + setNewParticleIDs(ion_tile, np_ion, num_added); + } +} diff --git a/Source/Particles/Collision/CMakeLists.txt b/Source/Particles/Collision/CMakeLists.txt index 7a4bc4c14..b6c4e928f 100644 --- a/Source/Particles/Collision/CMakeLists.txt +++ b/Source/Particles/Collision/CMakeLists.txt @@ -2,5 +2,7 @@ target_sources(WarpX PRIVATE CollisionHandler.cpp CollisionBase.cpp + BackgroundMCCCollision.cpp + MCCProcess.cpp ) diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp index 600bbff29..40889ccb0 100644 --- a/Source/Particles/Collision/CollisionHandler.cpp +++ b/Source/Particles/Collision/CollisionHandler.cpp @@ -8,6 +8,7 @@ #include "BinaryCollision.H" #include "PairWiseCoulombCollisionFunc.H" +#include "BackgroundMCCCollision.H" #include <AMReX_ParmParse.H> @@ -37,6 +38,12 @@ CollisionHandler::CollisionHandler() allcollisions[i] = std::make_unique<BinaryCollision<PairWiseCoulombCollisionFunc>>(collision_names[i]); } + else if (type == "background_mcc") { + allcollisions[i] = std::make_unique<BackgroundMCCCollision>(collision_names[i]); + } + else{ + amrex::Abort("Unknown collision type."); + } } diff --git a/Source/Particles/Collision/MCCProcess.H b/Source/Particles/Collision/MCCProcess.H new file mode 100644 index 000000000..c4038e496 --- /dev/null +++ b/Source/Particles/Collision/MCCProcess.H @@ -0,0 +1,116 @@ +/* Copyright 2021 Modern Electron + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PARTICLES_COLLISION_MCCPROCESS_H_ +#define WARPX_PARTICLES_COLLISION_MCCPROCESS_H_ + +#include <AMReX_Vector.H> +#include <AMReX_RandomEngine.H> +#include <AMReX_GpuContainers.H> + +enum class MCCProcessType { + INVALID, + ELASTIC, + BACK, + CHARGE_EXCHANGE, + EXCITATION, + IONIZATION, +}; + +class MCCProcess +{ +public: + template <typename T> + using VectorType = amrex::Gpu::ManagedVector<T>; + + MCCProcess ( + const std::string& scattering_process, + const std::string& cross_section_file, + const amrex::Real energy + ); + + template <typename InputVector> + MCCProcess ( + const std::string& scattering_process, + const InputVector&& energies, + const InputVector&& sigmas, + const amrex::Real energy + ); + + void* operator new ( std::size_t count ); + void operator delete (void* ptr); + + /** Get the collision cross-section using a simple linear interpolator. If + * the energy value is lower (higher) than the given energy range the + * first (last) cross-section value is used. + * + * @param E_coll collision energy in eV + * + */ + AMREX_GPU_HOST_DEVICE + amrex::Real getCrossSection (amrex::Real E_coll) const + { + if (E_coll < m_energy_lo) { + return m_sigma_lo; + } else if (E_coll > m_energy_hi) { + return m_sigma_hi; + } else { + using std::floor; + using std::ceil; + // calculate index of bounding energy pairs + amrex::Real temp = (E_coll - m_energy_lo) / m_dE; + int idx_1 = floor(temp); + int idx_2 = ceil(temp); + + // linearly interpolate to the given energy value + temp -= idx_1; + return ( + m_sigmas_data[idx_1] + (m_sigmas_data[idx_2] - m_sigmas_data[idx_1]) * temp + ); + } + } + + /** Read the given cross-section data file to memory. + * + * @param cross_section_file the path to the file containing the cross- + * section data + * @param energies vector storing energy values in eV + * @param sigmas vector storing cross-section values + * + */ + static + void readCrossSectionFile ( + const std::string cross_section_file, + VectorType<amrex::Real>& energies, + VectorType<amrex::Real>& sigmas + ); + + static + void sanityCheckEnergyGrid ( + const VectorType<amrex::Real>& energies, + amrex::Real dE + ); + + MCCProcessType m_type; + amrex::Real m_energy_penalty; + +private: + + static + MCCProcessType parseProcessType(const std::string& process); + + void init(); + + VectorType<amrex::Real> m_energies; + amrex::Real* m_energies_data; + VectorType<amrex::Real> m_sigmas; + amrex::Real* m_sigmas_data; + + int m_grid_size; + amrex::Real m_energy_lo, m_energy_hi, m_sigma_lo, m_sigma_hi, m_dE; +}; + +#endif // WARPX_PARTICLES_COLLISION_MCCPROCESS_H_ diff --git a/Source/Particles/Collision/MCCProcess.cpp b/Source/Particles/Collision/MCCProcess.cpp new file mode 100644 index 000000000..75778ce2f --- /dev/null +++ b/Source/Particles/Collision/MCCProcess.cpp @@ -0,0 +1,129 @@ +/* Copyright 2021 Modern Electron + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#include "MCCProcess.H" +#include "WarpX.H" + +MCCProcess::MCCProcess ( + const std::string& scattering_process, + const std::string& cross_section_file, + const amrex::Real energy ) + : m_type(parseProcessType(scattering_process)) + , m_energy_penalty(energy) +{ + amrex::Print() << "Reading file " << cross_section_file << " for " + << scattering_process << " scattering cross-sections.\n"; + + // read the cross-section data file into memory + readCrossSectionFile(cross_section_file, m_energies, m_sigmas); + + init(); +} + +template <typename InputVector> +MCCProcess::MCCProcess ( + const std::string& scattering_process, + const InputVector&& energies, + const InputVector&& sigmas, + const amrex::Real energy ) + : m_type(parseProcessType(scattering_process)) + , m_energy_penalty(energy) +{ + using std::begin; + using std::end; + m_energies.insert(m_energies.begin(), begin(energies), end(energies)); + m_sigmas .insert(m_sigmas.begin(), begin(sigmas), end(sigmas)); + + init(); +} + +void* +MCCProcess::operator new ( std::size_t count ) +{ + return amrex::The_Managed_Arena()->alloc(count); +} + +void +MCCProcess::operator delete ( void* ptr ) +{ + amrex::The_Managed_Arena()->free(ptr); +} + +void +MCCProcess::init() +{ + m_energies_data = m_energies.data(); + m_sigmas_data = m_sigmas.data(); + + // save energy grid parameters for easy use + m_grid_size = m_energies.size(); + m_energy_lo = m_energies[0]; + m_energy_hi = m_energies[m_grid_size-1]; + m_sigma_lo = m_sigmas[0]; + m_sigma_hi = m_sigmas[m_grid_size-1]; + m_dE = (m_energy_hi - m_energy_lo)/(m_grid_size - 1.); + + // sanity check cross-section energy grid + sanityCheckEnergyGrid(m_energies, m_dE); + + // check that the cross-section is 0 at the energy cost if the energy + // cost is > 0 - this is to prevent the possibility of negative left + // over energy after a collision event + if (m_energy_penalty > 0) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (getCrossSection(m_energy_penalty) == 0), + "Cross-section > 0 at energy cost for collision." + ); + } +} + +MCCProcessType +MCCProcess::parseProcessType(const std::string& scattering_process) +{ + if (scattering_process == "elastic") { + return MCCProcessType::ELASTIC; + } else if (scattering_process == "back") { + return MCCProcessType::BACK; + } else if (scattering_process == "charge_exchange") { + return MCCProcessType::CHARGE_EXCHANGE; + } else if (scattering_process == "ionization") { + return MCCProcessType::IONIZATION; + } else if (scattering_process.find("excitation") != std::string::npos) { + return MCCProcessType::EXCITATION; + } else { + return MCCProcessType::INVALID; + } +} + +void +MCCProcess::readCrossSectionFile ( + const std::string cross_section_file, + MCCProcess::VectorType<amrex::Real>& energies, + MCCProcess::VectorType<amrex::Real>& sigmas ) +{ + std::ifstream infile(cross_section_file); + double energy, sigma; + while (infile >> energy >> sigma) { + energies.push_back(energy); + sigmas.push_back(sigma); + } +} + +void +MCCProcess::sanityCheckEnergyGrid ( + const VectorType<amrex::Real>& energies, + amrex::Real dE + ) +{ + // confirm that the input data for the cross-section was provided with + // equal energy steps, otherwise the linear interpolation will fail + for (unsigned i = 1; i < energies.size(); i++) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (std::abs(energies[i] - energies[i-1] - dE) < dE / 100.0), + "Energy grid not evenly spaced." + ); + } +} diff --git a/Source/Particles/Collision/MCCScattering.H b/Source/Particles/Collision/MCCScattering.H new file mode 100644 index 000000000..52c321ff7 --- /dev/null +++ b/Source/Particles/Collision/MCCScattering.H @@ -0,0 +1,307 @@ +/* Copyright 2021 Modern Electron + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef WARPX_PARTICLES_COLLISION_MCC_SCATTERING_H_ +#define WARPX_PARTICLES_COLLISION_MCC_SCATTERING_H_ + +#include "MCCProcess.H" +#include "Utils/WarpXConst.H" +#include <AMReX_Random.H> +#include <AMReX_REAL.H> + +/** @file + * + * This file contains the implementation of the scattering processes available + * in the MCC handling. + */ + +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void getRandomVector ( amrex::Real& x, amrex::Real& y, amrex::Real& z, + amrex::RandomEngine const& engine ) +{ + using std::sqrt; + using std::cos; + using std::sin; + using namespace amrex::literals; + + // generate random unit vector in 3 dimensions + // https://mathworld.wolfram.com/SpherePointPicking.html + amrex::Real const theta = amrex::Random(engine) * 2.0_rt * MathConst::pi; + z = 2.0_rt * amrex::Random(engine) - 1.0_rt; + amrex::Real const xy = sqrt(1_rt - z*z); + x = xy * cos(theta); + y = xy * sin(theta); +} + + +/** \brief Function to perform scattering of a particle that results in a + * random velocity vector with given magnitude. This is used in excitation and + * ionization events. + * + * @param[in/out] ux, uy, uz colliding particle's velocity + * @param[in] vp velocity magnitude of the colliding particle after collision. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void RandomizeVelocity ( amrex::ParticleReal& ux, amrex::ParticleReal& uy, + amrex::ParticleReal& uz, + const amrex::ParticleReal vp, + amrex::RandomEngine const& engine ) +{ + amrex::Real x, y, z; + // generate random unit vector for the new velocity direction + getRandomVector(x, y, z, engine); + + // scale new vector to have the desired magnitude + ux = x * vp; + uy = y * vp; + uz = z * vp; +} + + +/** \brief Function to perform elastic scattering of a particle in the lab + * frame. The particle velocities transformed to the COM frame where a hard + * sphere collision occurs. The resulting particle velocities are transformed + * back to the lab frame and the input particle's velocity is updated. + * + * @param[in/out] ux, uy, uz colliding particle's velocity + * @param[in] uCOM_x, uCOM_y, uCOM_z velocity of the center of momentum frame. + * @param[in] engine random number generator. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ElasticScattering ( amrex::ParticleReal& ux, amrex::ParticleReal& uy, + amrex::ParticleReal& uz, amrex::ParticleReal uCOM_x, + amrex::ParticleReal uCOM_y, amrex::ParticleReal uCOM_z, + amrex::RandomEngine const& engine ) +{ + using std::sqrt; + + // transform to center of momentum frame + ux -= uCOM_x; + uy -= uCOM_y; + uz -= uCOM_z; + + // istropically scatter the particle + amrex::Real const mag = sqrt(ux*ux + uy*uy + uz*uz); + RandomizeVelocity(ux, uy, uz, mag, engine); + + // transform back to lab frame + ux += uCOM_x; + uy += uCOM_y; + uz += uCOM_z; +} + + +/** \brief Function to perform back scattering of a particle in the lab + * frame. The particle velocity is transformed to the COM frame where it is + * reversed. The resulting particle velocities are then transformed back to the + * lab frame and the input particle's velocity is updated. + * + * @param[in/out] ux, uy, uz colliding particle's velocity + * @param[in] uCOM_x, uCOM_y, uCOM_z velocity of the center of momentum frame. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void BackScattering ( amrex::ParticleReal& ux, amrex::ParticleReal& uy, + amrex::ParticleReal& uz, + const amrex::ParticleReal uCOM_x, + const amrex::ParticleReal uCOM_y, + const amrex::ParticleReal uCOM_z) +{ + // transform to COM frame, reverse particle velocity and transform back + using namespace amrex::literals; + ux = -1.0_prt * ux + 2.0_prt * uCOM_x; + uy = -1.0_prt * uy + 2.0_prt * uCOM_y; + uz = -1.0_prt * uz + 2.0_prt * uCOM_z; +} + + +/** \brief Function to perform charge exchange of an ion with a neutral + * particle. + * + * @param[in/out] ux, uy, uz colliding particle's velocity + * @param[in] ua_x, ua_y, ua_z velocity of the neutral particle. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ChargeExchange ( amrex::ParticleReal& ux, amrex::ParticleReal& uy, + amrex::ParticleReal& uz, + const amrex::ParticleReal ua_x, + const amrex::ParticleReal ua_y, + const amrex::ParticleReal ua_z) +{ + // swap ion velocity for neutral velocity + ux = ua_x; + uy = ua_y; + uz = ua_z; +} + + +/** + * \brief Filter functor for impact ionization + */ +class ImpactIonizationFilterFunc +{ +public: + + /** + * \brief Constructor of the ImpactIonizationFilterFunc functor. + * + * This function sample a random number and compares it to the total + * collision probability to see if the given particle should be considered + * for an ionization event. If the particle passes this stage the collision + * cross-section is calculated given its energy and another random number + * is used to determine whether it actually collides. + * + * @param[in] mcc_process an MCCProcess object associated with the ionization + * @param[in] mass colliding particle's mass (could also assume electron) + * @param[in] total_collision_prob total probability for a collision to occur + * @param[in] nu_max_norm normalized maximum collision frequency, normalized + * by the neutral density. + */ + ImpactIonizationFilterFunc( + MCCProcess mcc_process, + amrex::Real const mass, + amrex::Real const total_collision_prob, + amrex::Real const nu_max_norm + ) : m_mcc_process(mcc_process), m_mass(mass), + m_total_collision_prob(total_collision_prob), + m_nu_max_norm(nu_max_norm){ } + + /** + * \brief Functor call. This method determines if a given (electron) particle + * should undergo an ionization collision. + * + * @param[in] ptd particle tile data + * @param[in] i particle index + * @return true if a collision occurs, false otherwise + */ + template <typename PData> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + bool operator() ( + const PData& ptd, int const i, amrex::RandomEngine const& engine + ) const noexcept + { + using namespace amrex; + using std::sqrt; + + // determine if this particle should collide + if (Random(engine) > m_total_collision_prob) return false; + + // get the particle velocity + const ParticleReal ux = ptd.m_rdata[PIdx::ux][i]; + const ParticleReal uy = ptd.m_rdata[PIdx::uy][i]; + const ParticleReal uz = ptd.m_rdata[PIdx::uz][i]; + + // calculate kinetic energy + const ParticleReal v_coll2 = ux*ux + uy*uy + uz*uz; + const ParticleReal E_coll = 0.5_prt * m_mass * v_coll2 / PhysConst::q_e; + const ParticleReal v_coll = sqrt(v_coll2); + + // get collision cross-section + const Real sigma_E = m_mcc_process.getCrossSection(E_coll); + + // calculate normalized collision frequency + const Real nu_i = sigma_E * v_coll / m_nu_max_norm; + + // check if this collision should be performed + return (Random(engine) <= nu_i); + } + +private: + MCCProcess m_mcc_process; + amrex::Real m_mass; + amrex::Real m_total_collision_prob = 0; + amrex::Real m_nu_max_norm; +}; + + +/** + * \brief Transform functor for impact ionization + */ +class ImpactIonizationTransformFunc +{ +public: + + /** + * \brief Constructor of the ImpactIonizationTransformFunc functor. + * + * The transform is responsible for appropriately decreasing the kinetic + * energy of the colliding particle and assigning appropriate velocities + * to the two newly created particles. To this end the energy cost of + * ionization is passed to the constructor as well as the mass of the + * colliding species and the standard deviation of the ion velocity + * (normalized temperature). + * + * @param[in] energy_cost energy cost of ionization + * @param[in] mass1 mass of the colliding species + */ + ImpactIonizationTransformFunc( + amrex::Real energy_cost, amrex::Real mass1, amrex::Real ion_vel_std + ) : m_energy_cost(energy_cost), m_mass1(mass1), + m_ion_vel_std(ion_vel_std) { } + + /** + * \brief Functor call. It determines the properties of the generated pair + * and decreases the kinetic energy of the colliding particle. Inputs are + * standard from the FilterCopyTransfrom::filterCopyTransformParticles + * function. + * + * @param[in,out] dst1 target species 1 (electrons) + * @param[in,out] dst2 target species 2 (ions) + * @param[in] src source species (electrons) + * @param[in] i_src particle index of the source species + * @param[in] i_dst1 particle index of target species 1 + * @param[in] i_dst2 particle index of target species 2 + * @param[in] engine random number generator engine + */ + template <typename DstData, typename SrcData> + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void operator() (DstData& dst1, DstData& dst2, SrcData& src, + int const i_src, int const i_dst1, int const i_dst2, + amrex::RandomEngine const& engine) const noexcept + { + using namespace amrex; + using std::sqrt; + + // get references to the original particle's velocity + auto& ux = src.m_rdata[PIdx::ux][i_src]; + auto& uy = src.m_rdata[PIdx::uy][i_src]; + auto& uz = src.m_rdata[PIdx::uz][i_src]; + + // get references to the new particles' velocities + auto& e_ux = dst1.m_rdata[PIdx::ux][i_dst1]; + auto& e_uy = dst1.m_rdata[PIdx::uy][i_dst1]; + auto& e_uz = dst1.m_rdata[PIdx::uz][i_dst1]; + auto& i_ux = dst2.m_rdata[PIdx::ux][i_dst2]; + auto& i_uy = dst2.m_rdata[PIdx::uy][i_dst2]; + auto& i_uz = dst2.m_rdata[PIdx::uz][i_dst2]; + + // calculate kinetic energy + const ParticleReal v_coll2 = ux*ux + uy*uy + uz*uz; + const ParticleReal E_coll = 0.5_prt * m_mass1 * v_coll2 / PhysConst::q_e; + + // get the left over energy + amrex::Real E_remaining = E_coll - m_energy_cost; + + // each electron gets half the energy (could change this later) + amrex::Real vp = sqrt( + 2.0_prt / m_mass1 * PhysConst::q_e * E_remaining / 2.0_prt + ); + + // isotropically scatter electrons + RandomizeVelocity(ux, uy, uz, vp, engine); + RandomizeVelocity(e_ux, e_uy, e_uz, vp, engine); + + // get velocities for the ion from a Maxwellian distribution + i_ux = m_ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine); + i_uy = m_ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine); + i_uz = m_ion_vel_std * RandomNormal(0_prt, 1.0_prt, engine); + } + +private: + amrex::Real m_energy_cost; + amrex::Real m_mass1; + amrex::Real m_ion_vel_std; +}; +#endif // WARPX_PARTICLES_COLLISION_MCC_SCATTERING_H_ diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package index 4c7445cb6..4b01667d1 100644 --- a/Source/Particles/Collision/Make.package +++ b/Source/Particles/Collision/Make.package @@ -1,4 +1,6 @@ CEXE_sources += CollisionHandler.cpp CEXE_sources += CollisionBase.cpp +CEXE_sources += BackgroundMCCCollision.cpp +CEXE_sources += MCCProcess.cpp VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision |