QuantumSyncEngineWrapper.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_quantum_sync_engine_wrapper_h_
#define WARPX_quantum_sync_engine_wrapper_h_
 
#include "BreitWheelerEngineWrapper_fwd.H"
 
#include "QedChiFunctions.H"
#include "QedWrapperCommons.H"
#include "Utils/WarpXConst.H"
 
#include <AMReX_Extension.H>
#include <AMReX_GpuQualifiers.H>
#include <AMReX_REAL.H>
#include <AMReX_Random.H>
 
#include <picsar_qed/containers/picsar_array.hpp>
#include <picsar_qed/math/cmath_overloads.hpp>
#include <picsar_qed/math/vec_functions.hpp>
#include <picsar_qed/physics/phys_constants.h>
#include <picsar_qed/physics/quantum_sync/quantum_sync_engine_core.hpp>
#include <picsar_qed/physics/quantum_sync/quantum_sync_engine_tables.hpp>
#include <picsar_qed/physics/unit_conversion.hpp>
 
#include <cmath>
#include <vector>
 
namespace amrex { struct RandomEngine; }
 
// Aliases =============================
using QS_dndt_table_params =
    picsar::multi_physics::phys::quantum_sync::
    dndt_lookup_table_params<amrex::ParticleReal>;
 
using QS_dndt_table =
    picsar::multi_physics::phys::quantum_sync::
    dndt_lookup_table<
    amrex::ParticleReal,
    PicsarQedVector<amrex::ParticleReal>>;
 
using QS_dndt_table_view = QS_dndt_table::view_type;
 
using QS_phot_em_table_params =
    picsar::multi_physics::phys::quantum_sync::
    photon_emission_lookup_table_params<amrex::ParticleReal>;
 
using QS_phot_em_table =
    picsar::multi_physics::phys::quantum_sync::
    photon_emission_lookup_table<
    amrex::ParticleReal,
    PicsarQedVector<amrex::ParticleReal>>;
 
using QS_phot_em_table_view = QS_phot_em_table::view_type;
 
struct PicsarQuantumSyncCtrl
{
    QS_dndt_table_params dndt_params;
    QS_phot_em_table_params phot_em_params;
};
 
// Functors ==================================
 
// These functors allow using the core elementary functions of the library.
// They are generated by a factory class (QuantumSynchrotronEngine, see below).
// They can be included in GPU kernels.
 
class QuantumSynchrotronGetOpticalDepth
{
public:
    QuantumSynchrotronGetOpticalDepth () = default;
 
    AMREX_GPU_HOST_DEVICE
    AMREX_FORCE_INLINE
    amrex::ParticleReal operator() (amrex::RandomEngine const& engine) const noexcept
    {
        namespace pxr_qs = picsar::multi_physics::phys::quantum_sync;
 
        //A random number in [0,1) should be provided as an argument.
        return pxr_qs::get_optical_depth(amrex::Random(engine));
    }
};
//____________________________________________
 
class QuantumSynchrotronEvolveOpticalDepth
{
public:
 
    QuantumSynchrotronEvolveOpticalDepth () = default;
 
    QuantumSynchrotronEvolveOpticalDepth (
        const QS_dndt_table_view table_view,
        const amrex::ParticleReal qs_minimum_chi_part):
        m_table_view{table_view}, m_qs_minimum_chi_part{qs_minimum_chi_part}{}
 
    AMREX_GPU_DEVICE
    AMREX_FORCE_INLINE
    int operator()(
        const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz,
        const amrex::ParticleReal ex, const amrex::ParticleReal ey, const amrex::ParticleReal ez,
        const amrex::ParticleReal bx, const amrex::ParticleReal by, const amrex::ParticleReal bz,
        const amrex::Real dt, amrex::ParticleReal& opt_depth) const noexcept
    {
        using namespace amrex::literals;
        namespace pxr_p = picsar::multi_physics::phys;
        namespace pxr_qs = picsar::multi_physics::phys::quantum_sync;
 
        constexpr amrex::ParticleReal m_e = PhysConst::m_e;
        constexpr amrex::ParticleReal inv_c2 = 1._rt/(PhysConst::c*PhysConst::c);
        const amrex::ParticleReal gamma = std::sqrt(1._rt + (ux*ux + uy*uy + uz*uz)*inv_c2);
        const auto energy = gamma*m_e*PhysConst::c*PhysConst::c;
 
        const auto chi_part = QedUtils::chi_ele_pos(
            m_e*ux, m_e*uy, m_e*uz, ex, ey, ez, bx, by, bz);
 
        if (chi_part < m_qs_minimum_chi_part)
            return 0;
 
        const auto is_out = pxr_qs::evolve_optical_depth<
            amrex::ParticleReal,
            QS_dndt_table_view,
            pxr_p::unit_system::SI>(
                energy, chi_part, dt, opt_depth, m_table_view);
 
        return is_out;
    }
 
private:
    QS_dndt_table_view m_table_view;
    amrex::ParticleReal m_qs_minimum_chi_part;
};
 
class QuantumSynchrotronPhotonEmission
{
public:
 
    QuantumSynchrotronPhotonEmission () = default;
 
    QuantumSynchrotronPhotonEmission (
        const QS_phot_em_table_view table_view):
        m_table_view{table_view}{}
 
    AMREX_GPU_DEVICE
    AMREX_FORCE_INLINE
    bool operator()(
    amrex::ParticleReal& ux, amrex::ParticleReal& uy,
    amrex::ParticleReal& uz,
    const amrex::ParticleReal ex, const amrex::ParticleReal ey,
    const amrex::ParticleReal ez, const amrex::ParticleReal bx,
    const amrex::ParticleReal by, const amrex::ParticleReal bz,
    amrex::ParticleReal& g_ux, amrex::ParticleReal& g_uy,
    amrex::ParticleReal& g_uz, amrex::RandomEngine const& engine) const noexcept
    {
        using namespace amrex;
        namespace pxr_m = picsar::multi_physics::math;
        namespace pxr_p = picsar::multi_physics::phys;
        namespace pxr_qs = picsar::multi_physics::phys::quantum_sync;
 
        const auto rand_zero_one_minus_epsi = amrex::Random(engine);
 
        constexpr ParticleReal me = PhysConst::m_e;
        constexpr ParticleReal one_over_me = 1._prt/me;
 
        // Particle momentum is stored as gamma * velocity.
        // Convert to m * gamma * velocity
        auto px = ux*me;
        auto py = uy*me;
        auto pz = uz*me;
 
        const auto chi_particle = QedUtils::chi_ele_pos(
            px, py, pz, ex, ey, ez, bx, by, bz);
 
        auto momentum_particle = pxr_m::vec3<amrex::ParticleReal>{px, py, pz};
        auto momentum_photon = pxr_m::vec3<amrex::ParticleReal>();
 
        const auto is_out = pxr_qs::generate_photon_update_momentum<
            amrex::ParticleReal,
            QS_phot_em_table_view,
            pxr_p::unit_system::SI>(
                chi_particle, momentum_particle,
                rand_zero_one_minus_epsi,
                m_table_view,
                momentum_photon);
 
        ux = momentum_particle[0]*one_over_me;
        uy = momentum_particle[1]*one_over_me;
        uz = momentum_particle[2]*one_over_me;
        g_ux = momentum_photon[0]*one_over_me;
        g_uy = momentum_photon[1]*one_over_me;
        g_uz = momentum_photon[2]*one_over_me;
 
        return is_out;
    }
 
private:
    QS_phot_em_table_view m_table_view;
 
};
 
// Factory class =============================
 
class QuantumSynchrotronEngine
{
public:
    QuantumSynchrotronEngine () = default;
 
    QuantumSynchrotronGetOpticalDepth build_optical_depth_functor ();
 
    QuantumSynchrotronEvolveOpticalDepth build_evolve_functor ();
 
    QuantumSynchrotronPhotonEmission build_phot_em_functor ();
 
    bool are_lookup_tables_initialized () const;
 
    std::vector<char> export_lookup_tables_data () const;
 
    bool init_lookup_tables_from_raw_data (const std::vector<char>& raw_data,
        const amrex::ParticleReal qs_minimum_chi_part);
 
    void init_builtin_tables(const amrex::ParticleReal qs_minimum_chi_part);
 
    void compute_lookup_tables (PicsarQuantumSyncCtrl ctrl,
        const amrex::ParticleReal qs_minimum_chi_part);
 
    PicsarQuantumSyncCtrl get_default_ctrl() const;
 
    amrex::ParticleReal get_minimum_chi_part() const;
 
private:
    bool m_lookup_tables_initialized = false;
 
    //Variables to store the minimum chi parameters to enable
    //Quantum Synchrotron process
    amrex::ParticleReal m_qs_minimum_chi_part;
 
    QS_dndt_table m_dndt_table;
    QS_phot_em_table m_phot_em_table;
 
    void init_builtin_dndt_table();
    void init_builtin_phot_em_table();
};
 
//============================================
 
#endif //WARPX_quantum_sync_engine_wrapper_h_