ReferenceParticle.H

Go to the documentation of this file.

/* Copyright 2022-2023 The Regents of the University of California, through Lawrence
 *           Berkeley National Laboratory (subject to receipt of any required
 *           approvals from the U.S. Dept. of Energy). All rights reserved.
 *
 * This file is part of ImpactX.
 *
 * Authors: Axel Huebl
 * License: BSD-3-Clause-LBNL
 */
#ifndef IMPACTX_REFERENCE_PARTICLE_H
#define IMPACTX_REFERENCE_PARTICLE_H
 
#include <ablastr/constant.H>
 
#include <AMReX_BLassert.H>
#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_Math.H>
#include <AMReX_REAL.H>
#include <AMReX_SmallMatrix.H>
 
#include <cmath>
#include <stdexcept>
#include <string>
 
 
namespace impactx
{
    struct RefPart
    {
        amrex::ParticleReal s = 0.0;  
        amrex::ParticleReal x = 0.0;  
        amrex::ParticleReal y = 0.0;  
        amrex::ParticleReal z = 0.0;  
        amrex::ParticleReal t = 0.0;  
        amrex::ParticleReal px = 0.0; 
        amrex::ParticleReal py = 0.0; 
        amrex::ParticleReal pz = 0.0; 
        amrex::ParticleReal pt = 0.0; 
        amrex::ParticleReal mass = 0.0; 
        amrex::ParticleReal charge = 0.0; 
        amrex::ParticleReal gyromagnetic_anomaly = 0.0; 
 
        amrex::ParticleReal sedge = 0.0;  
        amrex::SmallMatrix<amrex::ParticleReal, 6, 6, amrex::Order::F, 1> map{}; 

        RefPart
        copy () const
        {
            return RefPart{*this};
        }

        void
        reset (bool keep_mass=false, bool keep_charge=false)
        {
            auto old_mass = mass;
            auto old_charge = charge;
 
            RefPart zero{};
            *this = zero;
 
            if (keep_mass) { mass = old_mass; }
            if (keep_charge) { charge = old_charge; }
        }

        RefPart &
        set_species (std::string const & species_name)
        {
            constexpr amrex::ParticleReal m_e = ablastr::constant::SI::m_e_v<amrex::ParticleReal>;
            constexpr amrex::ParticleReal m_p = ablastr::constant::SI::m_p_v<amrex::ParticleReal>;
            constexpr amrex::ParticleReal MeV_inv_c2 = ablastr::constant::SI::MeV_inv_c2_v<amrex::ParticleReal>;
 
            amrex::ParticleReal qe;
            amrex::ParticleReal massE;
            amrex::ParticleReal g_anomaly;
 
            // [known_species]
            if (species_name == "electron") {
                qe = -1.0;
                massE = m_e / MeV_inv_c2;
                g_anomaly = 0.00115965218062;
            } else if (species_name == "positron") {
                qe = 1.0;
                massE = m_e / MeV_inv_c2;
                g_anomaly = 0.00115965218062;
            } else if (species_name == "proton") {
                qe = 1.0;
                massE = m_p / MeV_inv_c2;
                g_anomaly = 1.7928473446;
            } else if (species_name == "Hminus") {
                qe = -1.0;
                massE = (m_p + 2.0 * m_e) / MeV_inv_c2;  // neglect ~14 eV/c^2 binding energy
                g_anomaly = 1.7928473446;
            }
            // [/known_species]
            else {
                throw std::runtime_error(
                    "Unknown species: '" + species_name +
                    "'. Known species: electron, positron, proton, Hminus. "
                    "For other species, set charge, mass, and gyromagnetic anomaly "
                    "individually via set_charge_qe(), set_mass_MeV(), and "
                    "set_gyromagnetic_anomaly()."
                );
            }
 
            set_charge_qe(qe);
            set_mass_MeV(massE);
            set_gyromagnetic_anomaly(g_anomaly);
            return *this;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        gamma () const
        {
            amrex::ParticleReal const ref_gamma = -pt;
            return ref_gamma;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        beta () const
        {
            using namespace amrex::literals;
            using amrex::Math::powi;
 
            amrex::ParticleReal const ref_gamma = -pt;
            amrex::ParticleReal const ref_beta = std::sqrt(1.0_prt - 1.0_prt / powi<2>(ref_gamma));
            return ref_beta;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        beta_gamma () const
        {
            using namespace amrex::literals;
            using amrex::Math::powi;
 
            amrex::ParticleReal const ref_gamma = -pt;
            amrex::ParticleReal const ref_betagamma = std::sqrt(powi<2>(ref_gamma) - 1.0_prt);
            return ref_betagamma;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        mass_MeV () const
        {
            using namespace amrex::literals;
 
            constexpr amrex::ParticleReal inv_MeV_invc2 = 1.0_prt / ablastr::constant::SI::MeV_inv_c2_v<amrex::ParticleReal>;
            return amrex::ParticleReal(mass * inv_MeV_invc2);
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        RefPart &
        set_mass_MeV (amrex::ParticleReal const massE)
        {
            using namespace amrex::literals;
            using amrex::Math::powi;
 
            AMREX_ASSERT_WITH_MESSAGE(massE != 0.0_prt,
                                      "set_mass_MeV: Mass cannot be zero!");
 
            mass = massE * ablastr::constant::SI::MeV_inv_c2_v<amrex::ParticleReal>;
 
            // re-scale pt and pz
            if (pt != 0.0_prt)
            {
                pt = -kin_energy_MeV() / massE - 1.0_prt;
                pz = std::sqrt(powi<2>(pt) - 1.0_prt);
            }
 
            return *this;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        kin_energy_MeV () const
        {
            using namespace amrex::literals;
 
            amrex::ParticleReal const ref_gamma = -pt;
            amrex::ParticleReal const ref_kin_energy = mass_MeV() * (ref_gamma - 1.0_prt);
            return ref_kin_energy;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        RefPart &
        set_kin_energy_MeV (amrex::ParticleReal const kin_energy)
        {
            using namespace amrex::literals;
            using amrex::Math::powi;
 
            AMREX_ASSERT_WITH_MESSAGE(mass != 0.0_prt,
                                      "set_kin_energy_MeV: Set mass first!");
 
            px = 0.0;
            py = 0.0;
            pt = -kin_energy / mass_MeV() - 1.0_prt;
            pz = std::sqrt(powi<2>(pt) - 1.0_prt);
 
            return *this;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        rigidity_Tm () const
        {
            using namespace amrex::literals;
            using amrex::Math::powi;
 
            amrex::ParticleReal const ref_gamma = -pt;
            amrex::ParticleReal const ref_betagamma = std::sqrt(powi<2>(ref_gamma) - 1.0_prt);
            amrex::ParticleReal const ref_rigidity = mass*ref_betagamma*(ablastr::constant::SI::c)/charge;
            return ref_rigidity;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        charge_qe () const
        {
            using namespace amrex::literals;
 
            constexpr amrex::ParticleReal inv_qe = 1.0_prt / ablastr::constant::SI::q_e;
            return amrex::ParticleReal(charge * inv_qe);
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        RefPart &
        set_charge_qe (amrex::ParticleReal const charge_qe)
        {
            using namespace amrex::literals;
 
            this->charge = charge_qe * ablastr::constant::SI::q_e;
 
            return *this;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        amrex::ParticleReal
        qm_ratio_SI () const
        {
            return charge / mass;
        }

        AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
        RefPart &
        set_gyromagnetic_anomaly (amrex::ParticleReal const gyromagnetic_anomaly_ref)
        {
            using namespace amrex::literals;
 
            this->gyromagnetic_anomaly = gyromagnetic_anomaly_ref;
 
            return *this;
        }
 
    };
 
} // namespace impactx
 
#endif // IMPACTX_REFERENCE_PARTICLE_H