diff options
Diffstat (limited to 'Source')
| -rw-r--r-- | Source/Make.WarpX | 17 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.H | 14 | ||||
| -rw-r--r-- | Source/Particles/MultiParticleContainer.cpp | 59 | ||||
| -rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 1 | ||||
| -rw-r--r-- | Source/QED/BreitWheelerEngineTableBuilder.H | 11 | ||||
| -rw-r--r-- | Source/QED/BreitWheelerEngineTableBuilder.cpp | 14 | ||||
| -rw-r--r-- | Source/QED/BreitWheelerEngineWrapper.H | 176 | ||||
| -rw-r--r-- | Source/QED/BreitWheelerEngineWrapper.cpp | 199 | ||||
| -rw-r--r-- | Source/QED/Make.package | 11 | ||||
| -rw-r--r-- | Source/QED/QedChiFunctions.H | 38 | ||||
| -rw-r--r-- | Source/QED/QedTableParserHelperFunctions.H | 44 | ||||
| -rw-r--r-- | Source/QED/QedWrapperCommons.H | 32 | ||||
| -rw-r--r-- | Source/QED/QuantumSyncEngineTableBuilder.H | 6 | ||||
| -rw-r--r-- | Source/QED/QuantumSyncEngineTableBuilder.cpp | 13 | ||||
| -rw-r--r-- | Source/QED/QuantumSyncEngineWrapper.H | 168 | ||||
| -rw-r--r-- | Source/QED/QuantumSyncEngineWrapper.cpp | 149 |
16 files changed, 930 insertions, 22 deletions
diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 54baecbf6..b9fe0a258 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -83,11 +83,6 @@ ifeq ($(USE_SENSEI_INSITU),TRUE) endif ifeq ($(QED),TRUE) - BOOST_ROOT ?= NOT_SET - ifneq ($(BOOST_ROOT),NOT_SET) - VPATH_LOCATIONS += $(BOOST_ROOT) - INCLUDE_LOCATIONS += $(BOOST_ROOT) - endif INCLUDE_LOCATIONS += $(PICSAR_HOME)/src/multi_physics/QED/src CXXFLAGS += -DWARPX_QED CFLAGS += -DWARPX_QED @@ -95,8 +90,20 @@ ifeq ($(QED),TRUE) F90FLAGS += -DWARPX_QED include $(WARPX_HOME)/Source/QED/Make.package USERSuffix := $(USERSuffix).QED + + ifeq ($(QED_TABLE_GEN),TRUE) + BOOST_ROOT ?= NOT_SET + ifneq ($(BOOST_ROOT),NOT_SET) + VPATH_LOCATIONS += $(BOOST_ROOT) + INCLUDE_LOCATIONS += $(BOOST_ROOT) + CXXFLAGS += -DWARPX_QED_TABLE_GEN + CFLAGS += -DWARPX_QED_TABLE_GEN + FFLAGS += -DWARPX_QED_TABLE_GEN + F90FLAGS += -DWARPX_QED_TABLE_GEN + endif endif +endif include $(PICSAR_HOME)/src/Make.package diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 58546a106..9f52dd0a5 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -3,6 +3,7 @@ #define WARPX_ParticleContainer_H_ #include <AMReX_Particles.H> +#include <AMReX_ParallelDescriptor.H> #include <WarpXParticleContainer.H> #include <PhysicalParticleContainer.H> #include <RigidInjectedParticleContainer.H> @@ -18,6 +19,8 @@ #include <map> #include <string> #include <algorithm> +#include <utility> +#include <tuple> // // MultiParticleContainer holds multiple (nspecies or npsecies+1 when @@ -222,6 +225,17 @@ protected: //Initialize QED engines and provides smart pointers //to species who need QED processes void InitQED (); + + bool someone_has_quantum_sync = false; + bool someone_has_breit_wheeler = false; + + void InitQuantumSync (); + void InitBreitWheeler (); + + std::tuple<bool, std::string, PicsarQuantumSynchrotronCtrl> + ParseQuantumSyncParams (); + std::tuple<bool, std::string, PicsarBreitWheelerCtrl> + ParseBreitWheelerParams (); #endif private: diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index c860d21f5..afeac1abd 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -5,7 +5,6 @@ #include <limits> #include <algorithm> #include <string> -#include <memory> using namespace amrex; @@ -740,11 +739,69 @@ void MultiParticleContainer::InitQED () if(pc->has_quantum_sync()){ pc->set_quantum_sync_engine_ptr (std::make_shared<QuantumSynchrotronEngine>(qs_engine)); + someone_has_quantum_sync = true; } if(pc->has_breit_wheeler()){ pc->set_breit_wheeler_engine_ptr (std::make_shared<BreitWheelerEngine>(bw_engine)); + someone_has_breit_wheeler = true; } } + + if(someone_has_quantum_sync) + InitQuantumSync(); + + if(someone_has_breit_wheeler) + InitBreitWheeler(); + +} + +void MultiParticleContainer::InitQuantumSync () +{ + bool is_custom; + PicsarQuantumSynchrotronCtrl ctrl; + std::string filename; + std::tie(is_custom, filename, ctrl) = ParseQuantumSyncParams(); + + // if(ParallelDescriptor::IOProcessor()){ + // qs_engine.compute_lookup_tables_default(); + // qs_engine.write_lookup_tables("qed_qs_lookup.bin"); + // } + // amrex::ParallelDescriptor::Barrier(); + // qs_engine.read_lookup_tables("qed_qs_lookup.bin"); +} + +void MultiParticleContainer::InitBreitWheeler () +{ + bool is_custom; + PicsarBreitWheelerCtrl ctrl; + std::string filename; + std::tie(is_custom, filename, ctrl) = ParseBreitWheelerParams(); + + // if(ParallelDescriptor::IOProcessor()){ + // bw_engine.compute_lookup_tables_default(); + // bw_engine.write_lookup_tables("qed_bw_lookup.bin"); + // } + // amrex::ParallelDescriptor::Barrier(); + // bw_engine.read_lookup_tables("qed_bw_lookup.bin"); + +} + +std::tuple<bool,std::string,PicsarQuantumSynchrotronCtrl> +MultiParticleContainer::ParseQuantumSyncParams () +{ + PicsarQuantumSynchrotronCtrl ctrl; + bool is_custom{false}; + + return std::make_tuple(is_custom, std::string(""), ctrl); +} + +std::tuple<bool,std::string,PicsarBreitWheelerCtrl> +MultiParticleContainer::ParseBreitWheelerParams () +{ + PicsarBreitWheelerCtrl ctrl; + bool is_custom{false}; + + return std::make_tuple(is_custom, std::string(""), ctrl); } #endif diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c953aa2d7..08493edcd 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -10,6 +10,7 @@ #ifdef WARPX_QED #include <QuantumSyncEngineWrapper.H> #include <BreitWheelerEngineWrapper.H> + #include <QedChiFunctions.H> #endif #include <map> diff --git a/Source/QED/BreitWheelerEngineTableBuilder.H b/Source/QED/BreitWheelerEngineTableBuilder.H new file mode 100644 index 000000000..0e3e45064 --- /dev/null +++ b/Source/QED/BreitWheelerEngineTableBuilder.H @@ -0,0 +1,11 @@ +#ifndef WARPX_breit_wheeler_engine_table_builder_h_ +#define WARPX_breit_wheeler_engine_table_builder_h_ + +#include "QedWrapperCommons.H" + +class BreitWheelerEngineTableBuilder{ + public: + +}; + +#endif //WARPX_breit_wheeler_engine_table_builder_h_ diff --git a/Source/QED/BreitWheelerEngineTableBuilder.cpp b/Source/QED/BreitWheelerEngineTableBuilder.cpp new file mode 100644 index 000000000..e85bfd523 --- /dev/null +++ b/Source/QED/BreitWheelerEngineTableBuilder.cpp @@ -0,0 +1,14 @@ +#include "BreitWheelerEngineTableBuilder.H" + + +//Include the full BW engine with table generation support +//(after some consistency tests) +#ifdef PXRMP_CORE_ONLY + #error The Table Builder is incompatible with PXRMP_CORE_ONLY +#endif + +#ifdef __PICSAR_MULTIPHYSICS_BREIT_WHEELER_ENGINE__ + #warning breit_wheeler_engine.hpp should not have been included before reaching this point. +#endif +#include "breit_wheeler_engine.hpp" +//_______________________________________________ diff --git a/Source/QED/BreitWheelerEngineWrapper.H b/Source/QED/BreitWheelerEngineWrapper.H index 7cdc66b44..d42150154 100644 --- a/Source/QED/BreitWheelerEngineWrapper.H +++ b/Source/QED/BreitWheelerEngineWrapper.H @@ -3,12 +3,43 @@ #include "QedWrapperCommons.H" +#include <AMReX_Array.H> +#include <AMReX_Gpu.H> + +#include <string> +#include <fstream> + //BW ENGINE from PICSAR +//Include only essential parts of the library in this file #define PXRMP_CORE_ONLY -#include "breit_wheeler_engine.hpp" +#include <breit_wheeler_engine.hpp> + +#ifdef WARPX_QED_TABLE_GEN + #include "BreitWheelerEngineTableBuilder.H" +#endif + +using PicsarBreitWheelerEngine = picsar::multi_physics:: + breit_wheeler_engine<amrex::Real, QedUtils::DummyStruct>; -using WarpXBreitWheelerWrapper = - picsar::multi_physics::breit_wheeler_engine<amrex::Real, DummyStruct>; +using PicsarBreitWheelerCtrl = + picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real>; + +// Struct to hold engine data ================ + +struct BreitWheelerEngineInnards +{ + // Control parameters + PicsarBreitWheelerCtrl ctrl; + + //Lookup table data + amrex::Gpu::ManagedVector<amrex::Real> TTfunc_coords; + amrex::Gpu::ManagedVector<amrex::Real> TTfunc_data; + + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_1; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_2; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_data; + //______ +}; // Functors ================================== @@ -29,12 +60,111 @@ public: AMREX_FORCE_INLINE amrex::Real operator() () const { - return WarpXBreitWheelerWrapper:: + return PicsarBreitWheelerEngine:: internal_get_optical_depth(amrex::Random()); } }; //____________________________________________ +// Evolution of the optical depth (returns true if +// an event occurs) +class BreitWheelerEvolveOpticalDepth +{ +public: + BreitWheelerEvolveOpticalDepth( + BreitWheelerEngineInnards* t_innards): + p_innards{t_innards}{}; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + bool operator()( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real dt, amrex::Real& opt_depth) const + { + bool has_event_happend{false}; + amrex::Real dummy_lambda{1.0}; + amrex::Real unused_event_time{0.0}; + + const auto table = picsar::multi_physics::lookup_1d<amrex::Real> + (p_innards->TTfunc_data.size(), + p_innards->TTfunc_coords.data(), + p_innards->TTfunc_data.data()); + + PicsarBreitWheelerEngine:: + internal_evolve_opt_depth_and_determine_event( + px, py, pz, + ex, ey, ez, + bx, by, bz, + dt, opt_depth, + has_event_happend, unused_event_time, + dummy_lambda, + table, + p_innards->ctrl); + + return has_event_happend; + } + +private: + BreitWheelerEngineInnards* p_innards; +}; + +// Generates an electron-positron pair via the Breit-Wheeler process +// (returns false if errors occur) +class BreitWheelerGeneratePairs +{ +public: + BreitWheelerGeneratePairs( + BreitWheelerEngineInnards* t_innards): + p_innards{t_innards}{}; + + template <size_t sampling> + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + bool operator()( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real weight, + amrex::Real* e_px, amrex::Real* e_py, amrex::Real* e_pz, + amrex::Real* p_px, amrex::Real* p_py, amrex::Real* p_pz, + amrex::Real* e_weight, amrex::Real* p_weight) const + { + amrex::Real dummy_lambda{1.0}; + picsar::multi_physics::picsar_array<amrex::Real, sampling> + rand_zero_one_minus_epsi; + for(auto& el : rand_zero_one_minus_epsi) + el = amrex::Random(); + + const auto table = picsar::multi_physics::lookup_2d<amrex::Real> + (p_innards->cum_distrib_coords_1.size(), + p_innards->cum_distrib_coords_1.data(), + p_innards->cum_distrib_coords_2.size(), + p_innards->cum_distrib_coords_2.data(), + p_innards->cum_distrib_data.data()); + + bool stat = PicsarBreitWheelerEngine:: + internal_generate_breit_wheeler_pairs( + px, py, pz, + ex, ey, ez, + bx, by, bz, + weight, sampling, + e_px, e_py, e_pz, + p_px, p_py, p_pz, + e_weight, p_weight, + dummy_lambda, + table, + p_innards->ctrl, + rand_zero_one_minus_epsi.data()); + + return stat; + } + +private: + BreitWheelerEngineInnards* p_innards; +}; + // Factory class ============================= /** @@ -49,6 +179,44 @@ public: * \brief Builds the functor to initialize the optical depth */ BreitWheelerGetOpticalDepth build_optical_depth_functor (); + + /* \brief Builds the functor to evolve the optical depth */ + BreitWheelerEvolveOpticalDepth build_evolve_functor (); + + /* \brief Builds the functor to generate the pairs */ + BreitWheelerGeneratePairs build_pair_functor (); + + /* \brief Checks if lookup tables are properly initialized */ + bool are_lookup_tables_initialized () const; + + /* \brief Reads lookup tables from 'file' on disk */ + bool init_lookup_tables_from_raw_data (const std::vector<char>& raw_data); + + /* \brief Writes lookup tables on disk in 'file' + * return false if it fails. */ + std::vector<char> export_lookup_tables_data () const; + + /* \brief Computes the Lookup tables using the default settings + * provided by the PICSAR library */ + void compute_lookup_tables_default (); + + /* \brief Computes the Lookup tables using user-defined settings */ + void compute_custom_lookup_tables (PicsarBreitWheelerCtrl ctrl); + +private: + bool lookup_tables_initialized = false; + + BreitWheelerEngineInnards innards; + +#ifdef WARPX_QED_TABLE_GEN + + BreitWheelerEngineTableBuilder table_builder; + + //Private function which actually computes the lookup tables + void computes_lookup_tables ( + PicsarBreitWheelerCtrl ctrl); + +#endif }; //============================================ diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp index 33e6d48d7..327b76411 100644 --- a/Source/QED/BreitWheelerEngineWrapper.cpp +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -1,16 +1,209 @@ #include "BreitWheelerEngineWrapper.H" + +#include "QedTableParserHelperFunctions.H" + +#include <utility> + +using namespace std; +using namespace QedUtils; +using namespace amrex; + //This file provides a wrapper aroud the breit_wheeler engine //provided by the PICSAR library -using namespace picsar::multi_physics; - // Factory class ============================= BreitWheelerEngine::BreitWheelerEngine (){} -BreitWheelerGetOpticalDepth BreitWheelerEngine::build_optical_depth_functor () +//Builds the functor to initialize the optical depth +BreitWheelerGetOpticalDepth +BreitWheelerEngine::build_optical_depth_functor () { return BreitWheelerGetOpticalDepth(); } +//Builds the functor to evolve the optical depth +BreitWheelerEvolveOpticalDepth +BreitWheelerEngine::build_evolve_functor () +{ + AMREX_ALWAYS_ASSERT(lookup_tables_initialized); + + return BreitWheelerEvolveOpticalDepth(&innards); +} + +//Builds the functor to generate the pairs +BreitWheelerGeneratePairs +BreitWheelerEngine::build_pair_functor () +{ + AMREX_ALWAYS_ASSERT(lookup_tables_initialized); + + return BreitWheelerGeneratePairs(&innards); +} + +bool BreitWheelerEngine::are_lookup_tables_initialized () const +{ + return lookup_tables_initialized; +} + +/* \brief Reads lookup tables from 'file' on disk */ +bool +BreitWheelerEngine::init_lookup_tables_from_raw_data ( + const vector<char>& raw_data) +{ + const char* p_data = raw_data.data(); + const char* const p_last = &raw_data.back(); + bool is_ok; + + //Header (control parameters) + tie(is_ok, innards.ctrl.chi_phot_min, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_phot_tdndt_min, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_tdndt_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_phot_tdndt_max, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_tdndt_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_phot_tdndt_how_many, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_tdndt_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_phot_tpair_min, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_tpair_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_phot_tpair_max, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_tpair_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_phot_tpair_how_many, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_phot_tpair_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, innards.ctrl.chi_frac_tpair_how_many, p_data) = + parse_raw_data<decltype(innards.ctrl.chi_frac_tpair_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + //___________________________ + + //Data + vector<Real> tndt_coords(innards.ctrl.chi_phot_tdndt_how_many); + vector<Real> tndt_data(innards.ctrl.chi_phot_tdndt_how_many); + vector<Real> cum_tab_coords1(innards.ctrl.chi_phot_tpair_how_many); + vector<Real> cum_tab_coords2(innards.ctrl.chi_frac_tpair_how_many); + vector<Real> cum_tab_data(innards.ctrl.chi_phot_tpair_how_many* + innards.ctrl.chi_frac_tpair_how_many); + + tie(is_ok, tndt_coords, p_data) = + parse_raw_data_vec<Real>( + p_data, tndt_coords.size(), p_last); + if(!is_ok) return false; + innards.TTfunc_coords.assign(tndt_coords.begin(), tndt_coords.end()); + + tie(is_ok, tndt_data, p_data) = + parse_raw_data_vec<Real>( + p_data, tndt_data.size(), p_last); + if(!is_ok) return false; + innards.TTfunc_data.assign(tndt_data.begin(), tndt_data.end()); + + tie(is_ok, cum_tab_coords1, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_coords1.size(), p_last); + if(!is_ok) return false; + innards.cum_distrib_coords_1.assign( + cum_tab_coords1.begin(), cum_tab_coords1.end()); + + tie(is_ok, cum_tab_coords2, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_coords2.size(), p_last); + if(!is_ok) return false; + innards.cum_distrib_coords_2.assign( + cum_tab_coords2.begin(), cum_tab_coords2.end()); + + tie(is_ok, cum_tab_data, p_data) = + parse_raw_data_vec<Real>( + p_data, cum_tab_data.size(), p_last); + if(!is_ok) return false; + innards.cum_distrib_data.assign( + cum_tab_data.begin(), cum_tab_data.end()); + + //___________________________ + lookup_tables_initialized = true; + + return true; +} + + +#ifdef WARPX_QED_TABLE_GEN +//Initializes the Lookup tables using the default settings +//provided by the library +void BreitWheelerEngine::compute_lookup_tables_default () +{ + //A control parameters structure + //with the default values provided by the library + WarpXBreitWheelerWrapperCtrl ctrl_default; + + computes_lookup_tables(ctrl_default); + + lookup_tables_initialized = true; +} + +// Computes the Lookup tables using user-defined settings +void BreitWheelerEngine::compute_custom_lookup_tables ( + WarpXBreitWheelerWrapperCtrl ctrl) +{ + computes_lookup_tables(ctrl); + + lookup_tables_initialized = true; +} + + +/* \brief Writes lookup tables on disk in 'file' + * return false if it fails. */ +std::vector<char> export_lookup_tables_data () const +{ + if(!lookup_tables_initialized) + return std::vector<char>; + + //TODO + return std::vector<char>; +} + +//Private function which actually computes the lookup tables +void BreitWheelerEngine::computes_lookup_tables ( + PicsarBreitWheelerCtrl ctrl) +{ + //Lambda is not actually used if S.I. units are enabled + PicsarBreitWheelerEngine bw_engine( + std::move(QedUtils::DummyStruct()), 1.0, ctrl); + + bw_engine.compute_dN_dt_lookup_table(); + bw_engine.compute_cumulative_pair_table(); + + auto bw_innards_picsar = bw_engine.export_innards(); + + //Copy data in a GPU-friendly data-structure + innards.ctrl = bw_innards_picsar.bw_ctrl; + innards.TTfunc_coords.assign(bw_innards_picsar.TTfunc_table_coords_ptr, + bw_innards_picsar.TTfunc_table_coords_ptr + + bw_innards_picsar.TTfunc_table_coords_how_many); + innards.TTfunc_data.assign(bw_innards_picsar.TTfunc_table_data_ptr, + bw_innards_picsar.TTfunc_table_data_ptr + + bw_innards_picsar.TTfunc_table_data_how_many); + //____ +} + +#endif + //============================================ diff --git a/Source/QED/Make.package b/Source/QED/Make.package index d4bad3f18..23e50b9d1 100644 --- a/Source/QED/Make.package +++ b/Source/QED/Make.package @@ -1,8 +1,19 @@ CEXE_headers += QedWrapperCommons.H +CEXE_headers += QedChiFunctions.H +CEXE_headers += QedTableParserHelperFunctions.H CEXE_headers += BreitWheelerEngineWrapper.H CEXE_headers += QuantumSyncsEngineWrapper.H CEXE_sources += BreitWheelerEngineWrapper.cpp CEXE_sources += QuantumSyncEngineWrapper.cpp +#Table generation is enabled only if QED_TABLE_GEN is +#set to true +ifeq ($(QED_TABLE_GEN),TRUE) + CEXE_headers += BreitWheelerEngineTableBuilder.H + CEXE_headers += QuantumSyncEngineTableBuilder.H + CEXE_sources += BreitWheelerEngineTableBuilder.cpp + CEXE_sources += QuantumSyncEngineTableBuilder.cpp +endif + INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/QedChiFunctions.H b/Source/QED/QedChiFunctions.H new file mode 100644 index 000000000..b150aae6b --- /dev/null +++ b/Source/QED/QedChiFunctions.H @@ -0,0 +1,38 @@ +#ifndef WARPX_amrex_qed_chi_functions_h_ +#define WARPX_amrex_qed_chi_functions_h_ + +#include "QedWrapperCommons.H" + +#define PXRMP_CORE_ONLY +#include <chi_functions.hpp> + +namespace QedUtils{ + //Thin wrappers around functions used to calculate chi parameters + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real chi_photon( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz) + { + const amrex::Real dummy_lambda = 1.0; + return picsar::multi_physics::chi_photon( + px, py, pz, ex, ey, ez, bx, by, bz, dummy_lambda); + } + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + amrex::Real chi_lepton( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz) + { + const amrex::Real dummy_lambda = 1.0; + return picsar::multi_physics::chi_lepton( + px, py, pz, ex, ey, ez, bx, by, bz, dummy_lambda); + } + //_________ +}; + +#endif //WARPX_amrex_qed_chi_functions_h_
\ No newline at end of file diff --git a/Source/QED/QedTableParserHelperFunctions.H b/Source/QED/QedTableParserHelperFunctions.H new file mode 100644 index 000000000..a4c73adc4 --- /dev/null +++ b/Source/QED/QedTableParserHelperFunctions.H @@ -0,0 +1,44 @@ +#ifndef WARPX_amrex_qed_table_parser_helper_functions_h_ +#define WARPX_amrex_qed_table_parser_helper_functions_h_ + +//This file contains helper functions to parse a char* array +//into a lookup table + +#include <vector> +#include <tuple> + +namespace QedUtils{ + template <class T> + std::tuple<bool, std::vector<T>, const char*>parse_raw_data_vec( + const char* p_data, size_t how_many, const char* const p_last) + { + std::vector<T> res; + if(p_data + sizeof(T)*how_many > p_last) + return std::make_tuple(false, res, nullptr); + + auto r_data = reinterpret_cast<const T*>(p_data); + + res.assign(r_data, r_data + how_many); + + p_data += sizeof(T)*how_many; + return std::make_tuple(true, res, p_data); + } + + template <class T> + std::tuple<bool, T, const char*> parse_raw_data( + const char* p_data, const char* const p_last) + { + T res; + if(p_data + sizeof(T) > p_last) + return std::make_tuple(false, res, nullptr); + + auto r_data = reinterpret_cast<const T*>(p_data); + + res = *r_data; + + p_data += sizeof(T); + return std::make_tuple(true, res, p_data); + } +}; + +#endif //WARPX_amrex_qed_table_parser_helper_functions_h_ diff --git a/Source/QED/QedWrapperCommons.H b/Source/QED/QedWrapperCommons.H index 821034c06..a972b8869 100644 --- a/Source/QED/QedWrapperCommons.H +++ b/Source/QED/QedWrapperCommons.H @@ -1,16 +1,42 @@ #ifndef WARPX_amrex_qed_wrapper_commons_h_ #define WARPX_amrex_qed_wrapper_commons_h_ -//Common definitions for the QED library wrappers +//Common definitions for the QED library wrappers and table builders #include <AMReX_AmrCore.H> +#include <AMReX_GpuQualifiers.H> //Sets the decorator for GPU #define PXRMP_GPU AMREX_GPU_DEVICE //Sets SI units in the library #define PXRMP_WITH_SI_UNITS -//An empty data type -struct DummyStruct{}; +namespace QedUtils{ + //An empty data type + struct DummyStruct{}; +}; + +//Control structures used by the engines +struct PicsarBreitWheelerCtrl +{ + _REAL chi_phot_min = + static_cast<_REAL>(__breit_wheeler_min_chi_phot); + + _REAL chi_phot_tdndt_min = + static_cast<_REAL>(__breit_wheeler_min_tdndt_chi_phot); + _REAL chi_phot_tdndt_max = + static_cast<_REAL>(__breit_wheeler_max_tdndt_chi_phot); + size_t chi_phot_tdndt_how_many = + __breit_wheeler_how_many_tdndt_chi_phot; + + _REAL chi_phot_tpair_min = + static_cast<_REAL>(__breit_wheeler_min_tpair_chi_phot); + _REAL chi_phot_tpair_max = + static_cast<_REAL>(__breit_wheeler_max_tpair_chi_phot); + size_t chi_phot_tpair_how_many = + __breit_wheeler_how_many_tpair_chi_phot; + size_t chi_frac_tpair_how_many = + __breit_wheeler_chi_frac_tpair_how_many; +}; #endif //WARPX_amrex_qed_wrapper_commons_h_ diff --git a/Source/QED/QuantumSyncEngineTableBuilder.H b/Source/QED/QuantumSyncEngineTableBuilder.H new file mode 100644 index 000000000..f555f9212 --- /dev/null +++ b/Source/QED/QuantumSyncEngineTableBuilder.H @@ -0,0 +1,6 @@ +#ifndef WARPX_quantum_sync_engine_table_builder_h_ +#define WARPX_quantum_sync_engine_table_builder_h_ + +class QuantumSynchrotronEngineTableBuilder{}; + +#endif //WARPX_quantum_sync_engine_table_builder_h_ diff --git a/Source/QED/QuantumSyncEngineTableBuilder.cpp b/Source/QED/QuantumSyncEngineTableBuilder.cpp new file mode 100644 index 000000000..68dbe5a36 --- /dev/null +++ b/Source/QED/QuantumSyncEngineTableBuilder.cpp @@ -0,0 +1,13 @@ +#include "QuantumSyncEngineTableBuilder.H" + +//Include the full QS engine with table generation support +//(after some consistency tests) +#ifdef PXRMP_CORE_ONLY + #error The Table Builder is incompatible with PXRMP_CORE_ONLY +#endif + +#ifdef __PICSAR_MULTIPHYSICS_BREIT_WHEELER_ENGINE__ + #warning quantum_sync_engine.hpp should not have been included before reaching this point. +#endif +#include "quantum_sync_engine.hpp" +//_______________________________________________ diff --git a/Source/QED/QuantumSyncEngineWrapper.H b/Source/QED/QuantumSyncEngineWrapper.H index a285dc8b6..5c61894a6 100644 --- a/Source/QED/QuantumSyncEngineWrapper.H +++ b/Source/QED/QuantumSyncEngineWrapper.H @@ -3,12 +3,40 @@ #include "QedWrapperCommons.H" +#include <AMReX_Array.H> +#include <AMReX_Gpu.H> + + +#include <string> +#include <fstream> + //QS ENGINE from PICSAR +//Include only essential parts of the library in this file #define PXRMP_CORE_ONLY -#include "quantum_sync_engine.hpp" +#include <quantum_sync_engine.hpp> + +using PicsarQuantumSynchrotronEngine = picsar::multi_physics:: + quantum_synchrotron_engine<amrex::Real, QedUtils::DummyStruct>; + +using PicsarQuantumSynchrotronCtrl = + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real>; + +// Struct to hold engine data ================ -using WarpXQuantumSynchrotronWrapper = - picsar::multi_physics::quantum_synchrotron_engine<amrex::Real, DummyStruct>; +struct QuantumSynchrotronEngineInnards +{ + // Control parameters + PicsarQuantumSynchrotronCtrl ctrl; + + //Lookup table data + amrex::Gpu::ManagedDeviceVector<amrex::Real> KKfunc_coords; + amrex::Gpu::ManagedDeviceVector<amrex::Real> KKfunc_data; + + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_1; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_coords_2; + amrex::Gpu::ManagedVector<amrex::Real> cum_distrib_data; + //______ +}; // Functors ================================== @@ -29,12 +57,110 @@ public: AMREX_FORCE_INLINE amrex::Real operator() () const { - return WarpXQuantumSynchrotronWrapper:: + return PicsarQuantumSynchrotronEngine:: internal_get_optical_depth(amrex::Random()); } }; //____________________________________________ +// Evolution of the optical depth (returns true if +// an event occurs) +class QuantumSynchrotronEvolveOpticalDepth +{ +public: + QuantumSynchrotronEvolveOpticalDepth( + QuantumSynchrotronEngineInnards* t_innards): + p_innards{t_innards}{}; + + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + bool operator()( + amrex::Real px, amrex::Real py, amrex::Real pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real dt, amrex::Real& opt_depth) const + { + bool has_event_happend{false}; + amrex::Real dummy_lambda{1.0}; + amrex::Real unused_event_time{0.0}; + + const auto table = picsar::multi_physics + ::lookup_1d<amrex::Real>( + p_innards->KKfunc_data.size(), + p_innards->KKfunc_coords.data(), + p_innards->KKfunc_data.data()); + + PicsarQuantumSynchrotronEngine:: + internal_evolve_opt_depth_and_determine_event( + px, py, pz, + ex, ey, ez, + bx, by, bz, + dt, opt_depth, + has_event_happend, unused_event_time, + dummy_lambda, + table, + p_innards->ctrl); + + return has_event_happend; + } + +private: + QuantumSynchrotronEngineInnards* p_innards; +}; + +// Generates a photon via the Quantum Synchrotron process +// and updates momentum accordingly (returns false if errors occur) +class QuantumSynchrotronGeneratePhotonAndUpdateMomentum +{ +public: + QuantumSynchrotronGeneratePhotonAndUpdateMomentum( + QuantumSynchrotronEngineInnards* t_innards): + p_innards{t_innards}{}; + + template <size_t sampling> + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + bool operator()( + amrex::Real* px, amrex::Real* py, amrex::Real* pz, + amrex::Real ex, amrex::Real ey, amrex::Real ez, + amrex::Real bx, amrex::Real by, amrex::Real bz, + amrex::Real weight, + amrex::Real* g_px, amrex::Real* g_py, amrex::Real* g_pz, + amrex::Real* g_weight) const + { + amrex::Real dummy_lambda{1.0}; + amrex::GpuArray<amrex::Real, sampling> + rand_zero_one_minus_epsi; + for(auto& el : rand_zero_one_minus_epsi) + el = amrex::Random(); + + const auto table = picsar::multi_physics::lookup_2d<amrex::Real> + (p_innards->cum_distrib_coords_1.size(), + p_innards->cum_distrib_coords_1.data(), + p_innards->cum_distrib_coords_2.size(), + p_innards->cum_distrib_coords_2.data(), + p_innards->cum_distrib_data.data()); + + bool stat = PicsarQuantumSynchrotronEngine:: + internal_generate_photons_and_update_momentum( + *px, *py, *pz, + ex, ey, ez, + bx, by, bz, + weight, sampling, + *g_px, *g_py, *g_pz, + *g_weight, + dummy_lambda, + table, + p_innards->ctrl, + rand_zero_one_minus_epsi.data()); + + return stat; + } + +private: + QuantumSynchrotronEngineInnards* p_innards; +}; + // Factory class ============================= /** @@ -49,6 +175,40 @@ public: * \brief Builds the functor to initialize the optical depth */ QuantumSynchrotronGetOpticalDepth build_optical_depth_functor (); + + /* \brief Builds the functor to evolve the optical depth */ + QuantumSynchrotronEvolveOpticalDepth build_evolve_functor (); + + + /* \brief Checks if lookup tables are properly initialized */ + bool are_lookup_tables_initialized () const; + + /* \brief Reads lookup tables from 'file' on disk */ + void read_lookup_tables (std::string file); + +private: + bool lookup_tables_initialized = false; + + QuantumSynchrotronEngineInnards innards; + +#ifdef WARPX_QED_TABLE_GEN +public: + /* \brief Computes the Lookup tables using the default settings + * provided by the PICSAR library */ + void compute_lookup_tables_default (); + + /* \brief Computes the Lookup tables using user-defined settings */ + void compute_custom_lookup_tables (WarpXQuantumSynchrotronWrapperCtrl ctrl); + /* \brief Writes lookup tables on disk in 'file' + * return false if it fails. */ + bool write_lookup_tables (std::string file) const; + +private: + //Private function which actually computes the lookup tables + void computes_lookup_tables ( + WarpXQuantumSynchrotronWrapperCtrl ctrl); + +#endif }; //============================================ diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp index b55149187..c362bded3 100644 --- a/Source/QED/QuantumSyncEngineWrapper.cpp +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -2,16 +2,161 @@ //This file provides a wrapper aroud the quantum_sync engine //provided by the PICSAR library -using namespace picsar::multi_physics; - // Factory class ============================= QuantumSynchrotronEngine::QuantumSynchrotronEngine (){} +//Builds the functor to evolve the optical depth QuantumSynchrotronGetOpticalDepth QuantumSynchrotronEngine::build_optical_depth_functor () { return QuantumSynchrotronGetOpticalDepth(); } +//Builds the functor to evolve the optical depth +QuantumSynchrotronEvolveOpticalDepth QuantumSynchrotronEngine::build_evolve_functor () +{ + AMREX_ALWAYS_ASSERT(lookup_tables_initialized); + + return QuantumSynchrotronEvolveOpticalDepth(&innards); +} + +bool QuantumSynchrotronEngine::are_lookup_tables_initialized () const +{ + return lookup_tables_initialized; +} + +/* \brief Reads lookup tables from 'file' on disk */ +void QuantumSynchrotronEngine::read_lookup_tables (std::string file) +{ + std::ifstream ifile(file, std::ios::in | std::ios::binary); + + //Header (control parameters) + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_min), + sizeof(innards.ctrl.chi_part_min)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_tdndt_min), + sizeof(innards.ctrl.chi_part_tdndt_min)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_tdndt_max), + sizeof(innards.ctrl.chi_part_tdndt_max)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_tdndt_how_many), + sizeof(innards.ctrl.chi_part_tdndt_how_many)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_tem_min), + sizeof(innards.ctrl.chi_part_tem_min)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_tem_max), + sizeof(innards.ctrl.chi_part_tem_max)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.chi_part_tem_how_many), + sizeof(innards.ctrl.chi_part_tem_how_many)); + ifile.read(reinterpret_cast<char*>(&innards.ctrl.prob_tem_how_many), + sizeof(innards.ctrl.prob_tem_how_many)); + //_______ + + //Data + size_t size_buf = sizeof(amrex::Real)*innards.ctrl.chi_part_tdndt_how_many; + char* data_buf = new char(size_buf); + ifile.read(data_buf, size_buf); + innards.KKfunc_coords.assign((amrex::Real*)data_buf, + (amrex::Real*)data_buf + size_buf); + ifile.read(data_buf, size_buf); + innards.KKfunc_data.assign((amrex::Real*)data_buf, + (amrex::Real*)data_buf + size_buf); + delete[] data_buf; + //_______ + + ifile.close(); + + lookup_tables_initialized = true; +} + +#ifdef WARPX_QED_TABLE_GEN + +//Initializes the Lookup tables using the default settings +//provided by the library +void QuantumSynchrotronEngine::compute_lookup_tables_default () +{ + //A control parameters structure + //with the default values provided by the library + WarpXQuantumSynchrotronWrapperCtrl ctrl_default; + + computes_lookup_tables(ctrl_default); + + lookup_tables_initialized = true; +} + +// Computes the Lookup tables using user-defined settings +void QuantumSynchrotronEngine::compute_custom_lookup_tables ( + WarpXQuantumSynchrotronWrapperCtrl ctrl) +{ + computes_lookup_tables(ctrl); + + lookup_tables_initialized = true; +} + + +/* \brief Writes lookup tables on disk in 'file' + * return false if it fails. */ +bool QuantumSynchrotronEngine::write_lookup_tables ( + std::string file) const +{ + if(!lookup_tables_initialized) + return false; + + std::ofstream of(file, std::ios::out | std::ios::binary); + + //Header (control parameters) + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_min), + sizeof(innards.ctrl.chi_part_min)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_tdndt_min), + sizeof(innards.ctrl.chi_part_tdndt_min)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_tdndt_max), + sizeof(innards.ctrl.chi_part_tdndt_max)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_tdndt_how_many), + sizeof(innards.ctrl.chi_part_tdndt_how_many)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_tem_min), + sizeof(innards.ctrl.chi_part_tem_max)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_tem_max), + sizeof(innards.ctrl.chi_part_tem_max)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.chi_part_tem_how_many), + sizeof(innards.ctrl.chi_part_tem_how_many)); + of.write(reinterpret_cast<const char*>(&innards.ctrl.prob_tem_how_many), + sizeof(innards.ctrl.prob_tem_how_many)); + //_______ + + //Data + of.write(reinterpret_cast<const char*>(innards.KKfunc_coords.dataPtr()), + sizeof(amrex::Real)*innards.KKfunc_coords.size()); + of.write(reinterpret_cast<const char*>(innards.KKfunc_data.dataPtr()), + sizeof(amrex::Real)*innards.KKfunc_data.size()); + // TODO: add other table + //_______ + + of.close(); + + return true; +} + + //Private function which actually computes the lookup tables +void QuantumSynchrotronEngine::computes_lookup_tables ( + WarpXQuantumSynchrotronWrapperCtrl ctrl) +{ + //Lambda is not actually used if S.I. units are enabled + WarpXQuantumSynchrotronWrapper qs_engine( + std::move(QedUtils::DummyStruct()), 1.0, ctrl); + + qs_engine.compute_dN_dt_lookup_table(); + //qs_engine.compute_cumulative_pair_table(); + + auto qs_innards_picsar = qs_engine.export_innards(); + + //Copy data in a GPU-friendly data-structure + innards.ctrl = qs_innards_picsar.qs_ctrl; + innards.KKfunc_coords.assign(qs_innards_picsar.KKfunc_table_coords_ptr, + qs_innards_picsar.KKfunc_table_coords_ptr + + qs_innards_picsar.KKfunc_table_coords_how_many); + innards.KKfunc_data.assign(qs_innards_picsar.KKfunc_table_data_ptr, + qs_innards_picsar.KKfunc_table_data_ptr + + qs_innards_picsar.KKfunc_table_data_how_many); + //____ +} +#endif + //============================================ |