QEDPairGeneration.H

Go to the documentation of this file.

/* Copyright 2019 Luca Fedeli
 *
 * This file is part of WarpX.
 *
 * License: BSD-3-Clause-LBNL
 */
 
#ifndef WARPX_QED_PAIR_GENERATION_H_
#define WARPX_QED_PAIR_GENERATION_H_
 
#include "Particles/Gather/FieldGather.H"
#include "Particles/Gather/GetExternalFields.H"
#include "Particles/Pusher/GetAndSetPosition.H"
#include "Particles/WarpXParticleContainer.H"
#include "QEDInternals/BreitWheelerEngineWrapper.H"
 
#include <AMReX_Array.H>
#include <AMReX_Array4.H>
#include <AMReX_Dim3.H>
#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_IndexType.H>
#include <AMReX_REAL.H>
 
#include <AMReX_BaseFwd.H>


class PairGenerationFilterFunc
{
public:

    PairGenerationFilterFunc(int const opt_depth_runtime_comp)
        : m_opt_depth_runtime_comp(opt_depth_runtime_comp)
        {}

    template <typename PData>
    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
    bool operator() (const PData& ptd, int const i, amrex::RandomEngine const&) const noexcept
    {
        using namespace amrex;
 
        const amrex::ParticleReal opt_depth =
            ptd.m_runtime_rdata[m_opt_depth_runtime_comp][i];
        return (opt_depth < 0.0_rt);
    }
 
private:
    int m_opt_depth_runtime_comp = 0; 
};

class PairGenerationTransformFunc
{
public:

    PairGenerationTransformFunc(BreitWheelerGeneratePairs generate_functor,
                                const WarpXParIter& a_pti, int lev, amrex::IntVect ngEB,
                                amrex::FArrayBox const& exfab,
                                amrex::FArrayBox const& eyfab,
                                amrex::FArrayBox const& ezfab,
                                amrex::FArrayBox const& bxfab,
                                amrex::FArrayBox const& byfab,
                                amrex::FArrayBox const& bzfab,
                                amrex::Vector<amrex::ParticleReal>& E_external_particle,
                                amrex::Vector<amrex::ParticleReal>& B_external_particle,
                                int a_offset = 0);

    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;
 
        // gather E and B
        amrex::ParticleReal xp, yp, zp;
        m_get_position(i_src, xp, yp, zp);
 
        amrex::ParticleReal ex = m_Ex_external_particle;
        amrex::ParticleReal ey = m_Ey_external_particle;
        amrex::ParticleReal ez = m_Ez_external_particle;
        amrex::ParticleReal bx = m_Bx_external_particle;
        amrex::ParticleReal by = m_By_external_particle;
        amrex::ParticleReal bz = m_Bz_external_particle;
 
        m_get_externalEB(i_src, ex, ey, ez, bx, by, bz);
 
        doGatherShapeN(xp, yp, zp, ex, ey, ez, bx, by, bz,
                       m_ex_arr, m_ey_arr, m_ez_arr, m_bx_arr, m_by_arr, m_bz_arr,
                       m_ex_type, m_ey_type, m_ez_type, m_bx_type, m_by_type, m_bz_type,
                       m_dinv, m_xyzmin, m_lo, m_n_rz_azimuthal_modes,
                       m_nox, m_galerkin_interpolation);
 
        //Despite the names of the variables, positrons and electrons
        //can be exchanged, since the physical process is completely
        //symmetric with respect to this exchange.
        const auto& ux = src.m_rdata[PIdx::ux][i_src];
        const auto& uy = src.m_rdata[PIdx::uy][i_src];
        const auto& uz = src.m_rdata[PIdx::uz][i_src];
        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& p_ux = dst2.m_rdata[PIdx::ux][i_dst2];
        auto& p_uy = dst2.m_rdata[PIdx::uy][i_dst2];
        auto& p_uz = dst2.m_rdata[PIdx::uz][i_dst2];
        m_generate_functor(
            ux, uy, uz,
            ex, ey, ez,
            bx, by, bz,
            e_ux, e_uy, e_uz,
            p_ux, p_uy, p_uz,
            engine);
 
        src.m_idcpu[i_src] = amrex::ParticleIdCpus::Invalid;  // destroy photon after pair generation
    }
 
private:
 
    BreitWheelerGeneratePairs
    m_generate_functor; 
 
    GetParticlePosition<PIdx> m_get_position;
    GetExternalEBField m_get_externalEB;
    amrex::ParticleReal m_Ex_external_particle;
    amrex::ParticleReal m_Ey_external_particle;
    amrex::ParticleReal m_Ez_external_particle;
    amrex::ParticleReal m_Bx_external_particle;
    amrex::ParticleReal m_By_external_particle;
    amrex::ParticleReal m_Bz_external_particle;
 
    amrex::Array4<const amrex::Real> m_ex_arr;
    amrex::Array4<const amrex::Real> m_ey_arr;
    amrex::Array4<const amrex::Real> m_ez_arr;
    amrex::Array4<const amrex::Real> m_bx_arr;
    amrex::Array4<const amrex::Real> m_by_arr;
    amrex::Array4<const amrex::Real> m_bz_arr;
 
    amrex::IndexType m_ex_type;
    amrex::IndexType m_ey_type;
    amrex::IndexType m_ez_type;
    amrex::IndexType m_bx_type;
    amrex::IndexType m_by_type;
    amrex::IndexType m_bz_type;
 
    amrex::XDim3 m_dinv;
    amrex::XDim3 m_xyzmin;
 
    bool m_galerkin_interpolation;
    int m_nox;
    int m_n_rz_azimuthal_modes;
 
    amrex::Dim3 m_lo;
};
 
#endif //WARPX_QED_PAIR_GENERATION_H_