MultiParticleContainer.H

Go to the documentation of this file.

/* Copyright 2019-2020 Andrew Myers, Ann Almgren, Axel Huebl
 * David Grote, Jean-Luc Vay, Junmin Gu
 * Luca Fedeli, Mathieu Lobet, Maxence Thevenet
 * Remi Lehe, Revathi Jambunathan, Weiqun Zhang
 * Yinjian Zhao
 *
 * This file is part of WarpX.
 *
 * License: BSD-3-Clause-LBNL
 */
#ifndef WARPX_ParticleContainer_H_
#define WARPX_ParticleContainer_H_
 
#include "MultiParticleContainer_fwd.H"
#include "FieldSolver/ImplicitSolvers/ImplicitOptions.H"
 
#include "Particles/Collision/CollisionHandler.H"
#ifdef WARPX_QED
#   include "Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper_fwd.H"
#   include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper_fwd.H"
#endif
#include "PhysicalParticleContainer.H"
#include "Utils/TextMsg.H"
#include "Utils/WarpXConst.H"
#include "WarpXParticleContainer.H"
#include "ParticleBoundaries.H"
#include "ExternalParticleFields.H"
 
#include <ablastr/fields/MultiFabRegister.H>
 
#include <AMReX_BLassert.H>
#include <AMReX_Box.H>
#include <AMReX_Config.H>
#include <AMReX_Geometry.H>
#include <AMReX_GpuControl.H>
#include <AMReX_INT.H>
#include <AMReX_MFIter.H>
#include <AMReX_REAL.H>
#include <AMReX_RealBox.H>
#include <AMReX_Vector.H>
 
#include <AMReX_BaseFwd.H>
#include <AMReX_AmrCoreFwd.H>
 
#include <algorithm>
#include <array>
#include <iosfwd>
#include <iterator>
#include <limits>
#include <memory>
#include <string>
#include <vector>

class MultiParticleContainer
{
 
public:
 
    explicit MultiParticleContainer (amrex::AmrCore* amr_core);
 
    ~MultiParticleContainer() = default;
 
    MultiParticleContainer (MultiParticleContainer const &)             = delete;
    MultiParticleContainer& operator= (MultiParticleContainer const & ) = delete;
    MultiParticleContainer(MultiParticleContainer&& )                   = default;
    MultiParticleContainer& operator=(MultiParticleContainer&& )        = default;
 
    [[nodiscard]] WarpXParticleContainer&
    GetParticleContainer (int index) const {return *allcontainers[index];}
 
    [[nodiscard]] WarpXParticleContainer*
    GetParticleContainerPtr (int index) const {return allcontainers[index].get();}
 
    [[nodiscard]] WarpXParticleContainer&
    GetParticleContainerFromName (const std::string& name) const;
 
    std::array<amrex::ParticleReal, 3> meanParticleVelocity(int index) {
        return allcontainers[index]->meanParticleVelocity();
    }
 
    amrex::ParticleReal maxParticleDtInv();
 
    void TransformMomentumToCurvilinear (bool forward);
 
    void AllocData ();
 
    void InitData ();
 
    void InitMultiPhysicsModules ();

    void Evolve (
        ablastr::fields::MultiFabRegister& fields,
        int lev,
        std::string const& current_fp_string,
        amrex::Real t,
        amrex::Real dt,
        SubcyclingHalf subcycling_half=SubcyclingHalf::None,
        bool skip_deposition=false,
        PositionPushType position_push_type=PositionPushType::Full,
        MomentumPushType momentum_push_type=MomentumPushType::Full,
        ImplicitOptions const * implicit_options = nullptr
    );

    void DepositMassMatrices (ablastr::fields::MultiFabRegister& fields,
                              int lev, amrex::Real dt);

    void PushX (amrex::Real dt);

    void PushP (int lev, amrex::Real dt,
                const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
                const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz,
                MomentumPushType momentum_push_type);

    std::unique_ptr<amrex::MultiFab> GetZeroChargeDensity(int lev);

    void
    DepositCharge (const ablastr::fields::MultiLevelScalarField& rho,
                   amrex::Real relative_time);

    void
    DepositCurrent (ablastr::fields::MultiLevelVectorField const& J,
                    amrex::Real dt, amrex::Real relative_time,
                    PushType push_type = PushType::Explicit);

    void
    DepositTemperatures (ablastr::fields::MultiFabRegister & fields, amrex::Real relative_time);

 
    std::unique_ptr<amrex::MultiFab> GetChargeDensity(int lev, bool local = false);
 
    std::unique_ptr<amrex::MultiFab> GetGlobalPlasmaFrequency (int lev);
    void GenerateGlobalDebyeLength ();
 
    /*
     * \brief Calculates the electron-ion collision frequency for the specified species
     * @param[in] electron_species_name name of the electon species
     */
    void CalculateNuei (std::string const & electron_species_name);
 
    void doFieldIonization (int lev,
                            const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez,
                            const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz);
 
    void doCollisions (int step, amrex::Real cur_time, amrex::Real dt);

    void doResampling (const amrex::Vector<amrex::Geometry>& geom, int timestep, bool verbose);
 
#ifdef WARPX_QED
    void doQEDSchwinger ();

    [[nodiscard]] amrex::Box ComputeSchwingerGlobalBox () const;
#endif
 
    void Restart (const std::string& dir);
 
    void PostRestart ();
 
    void ReadHeader (std::istream& is);
 
    void WriteHeader (std::ostream& os) const;
 
    void SortParticlesByBin (
        const amrex::IntVect& bin_size,
        bool sort_particles_for_deposition,
        const amrex::IntVect& sort_idx_type);
 
    void Redistribute ();
 
    void defineAllParticleTiles ();
 
    void deleteInvalidParticles ();
 
    void RedistributeLocal (const amrex::IntVect& max_cells_travelled);

    void ApplyBoundaryConditions ();

    [[nodiscard]] amrex::Vector<amrex::Long> GetZeroParticlesInGrid(int lev) const;
 
    [[nodiscard]] amrex::Vector<amrex::Long> NumberOfParticlesInGrid(int lev) const;
 
    void Increment (amrex::MultiFab& mf, int lev);
 
    void SetParticleBoxArray (int lev, amrex::BoxArray& new_ba);
    void SetParticleDistributionMap (int lev, amrex::DistributionMapping& new_dm);
 
    [[nodiscard]] int nSpecies () const {return static_cast<int>(species_names.size());}
    [[nodiscard]] int nLasers () const {return static_cast<int>(lasers_names.size());}
    [[nodiscard]] int nContainers () const {return static_cast<int>(allcontainers.size());}

    void SetDoBackTransformedParticles (bool do_back_transformed_particles);
    void SetDoBackTransformedParticles (const std::string& species_name, bool do_back_transformed_particles);

    [[nodiscard]] bool getDoBackTransformedParticles () const
    {
        return m_do_back_transformed_particles;
    }
 
    [[nodiscard]] int nSpeciesDepositOnMainGrid () const
    {
        bool const onMainGrid = true;
        auto const & v = m_deposit_on_main_grid;
        return static_cast<int>(std::count( v.begin(), v.end(), onMainGrid ));
    }
 
    [[nodiscard]] int nSpeciesGatherFromMainGrid() const
    {
        bool const fromMainGrid = true;
        auto const & v = m_gather_from_main_grid;
        return static_cast<int>(std::count( v.begin(), v.end(), fromMainGrid ));
    }
 
    // Inject particles during the simulation (for particles entering the
    // simulation domain after some iterations, due to flowing plasma and/or
    // moving window).
    void ContinuousInjection(const amrex::RealBox& injection_box) const;

    void UpdateAntennaPosition(amrex::Real dt) const;
 
    [[nodiscard]] int doContinuousInjection() const;
 
    // Inject particles from a surface during the simulation
    void ContinuousFluxInjection(amrex::Real t, amrex::Real dt) const;
 
    [[nodiscard]] std::vector<std::string> GetSpeciesNames() const { return species_names; }
 
    [[nodiscard]] std::vector<std::string> GetLasersNames() const { return lasers_names; }
 
    [[nodiscard]] std::vector<std::string> GetSpeciesAndLasersNames() const
    {
        std::vector<std::string> tmp = species_names;
        tmp.insert(tmp.end(), lasers_names.begin(), lasers_names.end());
        return tmp;
    }
 
    void ScrapeParticlesAtEB (ablastr::fields::MultiLevelScalarField const& distance_to_eb);
 
    std::string m_B_ext_particle_s = "none";
    std::string m_E_ext_particle_s = "none";
    // 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;
    std::unique_ptr<amrex::Parser> m_Bz_particle_parser;
    // Parser for E_external on the particle
    std::unique_ptr<amrex::Parser> m_Ex_particle_parser;
    std::unique_ptr<amrex::Parser> m_Ey_particle_parser;
    std::unique_ptr<amrex::Parser> m_Ez_particle_parser;
 
    ExternalParticleFields m_external_particle_fields_metadata;
 
    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
    void doQedEvents (int lev,
                      const amrex::MultiFab& Ex,
                      const amrex::MultiFab& Ey,
                      const amrex::MultiFab& Ez,
                      const amrex::MultiFab& Bx,
                      const amrex::MultiFab& By,
                      const amrex::MultiFab& Bz);
#endif
 
    [[nodiscard]] int getSpeciesID (const std::string& product_str) const;
 
    amrex::Vector<std::unique_ptr<WarpXParticleContainer>>::iterator begin() {return allcontainers.begin();}
    amrex::Vector<std::unique_ptr<WarpXParticleContainer>>::iterator end() {return allcontainers.end();}
 
protected:
 
#ifdef WARPX_QED
    void doQedBreitWheeler (int lev,
                            const amrex::MultiFab& Ex,
                            const amrex::MultiFab& Ey,
                            const amrex::MultiFab& Ez,
                            const amrex::MultiFab& Bx,
                            const amrex::MultiFab& By,
                            const amrex::MultiFab& Bz);

    void doQedQuantumSync (int lev,
                           const amrex::MultiFab& Ex,
                           const amrex::MultiFab& Ey,
                           const amrex::MultiFab& Ez,
                           const amrex::MultiFab& Bx,
                           const amrex::MultiFab& By,
                           const amrex::MultiFab& Bz);
#endif
 
    // Particle container types
    enum struct PCTypes {Physical, RigidInjected, Photon};
 
    std::vector<std::string> species_names;
 
    std::vector<std::string> lasers_names;
 
    std::unique_ptr<CollisionHandler> collisionhandler;

    std::vector<bool> m_deposit_on_main_grid;
    std::vector<bool> m_laser_deposit_on_main_grid;

    std::vector<bool> m_gather_from_main_grid;
 
    std::vector<PCTypes> species_types;
 
    template<typename ...Args>
    [[nodiscard]]
    amrex::MFItInfo getMFItInfo (const WarpXParticleContainer& pc_src,
                                 Args const&... pc_dsts) const noexcept
    {
        amrex::MFItInfo info;
 
        MFItInfoCheckTiling(pc_src, pc_dsts...);
 
        if (WarpXParticleContainer::do_tiling && amrex::Gpu::notInLaunchRegion()) {
            info.EnableTiling(WarpXParticleContainer::tile_size);
        }
 
#ifdef AMREX_USE_OMP
        info.SetDynamic(true);
#endif
 
        return info;
    }
 
 
#ifdef WARPX_QED
    // The QED engines
    std::shared_ptr<BreitWheelerEngine> m_shr_p_bw_engine;
    std::shared_ptr<QuantumSynchrotronEngine> m_shr_p_qs_engine;
    //_______________________________

    void InitQED ();
 
    //Variables to store how many species need a QED process
    int m_nspecies_quantum_sync = 0;
    int m_nspecies_breit_wheeler = 0;
    //________
 
    static constexpr auto m_default_quantum_sync_photon_creation_energy_threshold =
        static_cast<amrex::ParticleReal>(
        2.0 * PhysConst::m_e * PhysConst::c2 );  
 
 
    amrex::ParticleReal m_quantum_sync_photon_creation_energy_threshold =
        m_default_quantum_sync_photon_creation_energy_threshold; 

    [[nodiscard]] int NSpeciesQuantumSync() const { return  m_nspecies_quantum_sync;}

    [[nodiscard]] int NSpeciesBreitWheeler() const { return m_nspecies_breit_wheeler;}

    void InitQuantumSync ();

    void InitBreitWheeler ();

    void QuantumSyncGenerateTable();

    void BreitWheelerGenerateTable();

    bool m_do_qed_schwinger = false;
    std::string m_qed_schwinger_ele_product_name;
    std::string m_qed_schwinger_pos_product_name;
    int m_qed_schwinger_ele_product;
    int m_qed_schwinger_pos_product;
    amrex::Real m_qed_schwinger_y_size;
    int m_qed_schwinger_threshold_poisson_gaussian = 25;
    amrex::Real m_qed_schwinger_xmin = std::numeric_limits<amrex::Real>::lowest();
    amrex::Real m_qed_schwinger_xmax = std::numeric_limits<amrex::Real>::max();
    amrex::Real m_qed_schwinger_ymin = std::numeric_limits<amrex::Real>::lowest();
    amrex::Real m_qed_schwinger_ymax = std::numeric_limits<amrex::Real>::max();
    amrex::Real m_qed_schwinger_zmin = std::numeric_limits<amrex::Real>::lowest();
    amrex::Real m_qed_schwinger_zmax = std::numeric_limits<amrex::Real>::max();
 
#endif
 
private:
 
    // physical particles (+ laser)
    amrex::Vector<std::unique_ptr<WarpXParticleContainer>> allcontainers;
 
    void ReadParameters ();

    bool m_params_initialized = false;
 
    void mapSpeciesProduct ();
 
    bool m_do_back_transformed_particles = false;
 
    void MFItInfoCheckTiling(const WarpXParticleContainer& /*pc_src*/) const noexcept
    {}
 
    template<typename First, typename ...Args>
    void MFItInfoCheckTiling(const WarpXParticleContainer& pc_src,
                             First const& pc_dst, Args const&... others) const noexcept
    {
        if (WarpXParticleContainer::do_tiling && amrex::Gpu::notInLaunchRegion()) {
            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(pc_dst.do_tiling,
                "For particle creation processes, either all or none of the "
                "particle species must use tiling.");
        }
 
        MFItInfoCheckTiling(pc_src, others...);
    }

    void CheckIonizationProductSpecies();
 
#ifdef WARPX_QED
    void CheckQEDProductSpecies();
#endif
 
 
};
#endif /*WARPX_ParticleContainer_H_*/