diff options
Diffstat (limited to 'Source')
24 files changed, 1714 insertions, 71 deletions
diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index d55660b39..3b7481b8c 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -113,7 +113,7 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const } #ifdef WARPX_QED - if(pc->do_qed){ + if(pc->m_do_qed){ real_names.push_back("tau"); } #endif diff --git a/Source/Make.WarpX b/Source/Make.WarpX index 54baecbf6..b81fed147 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 -endif + ifeq ($(QED_TABLE_GEN),TRUE) + BOOST_ROOT ?= NOT_SET + ifneq ($(BOOST_ROOT),NOT_SET) + VPATH_LOCATIONS += $(BOOST_ROOT) + INCLUDE_LOCATIONS += $(BOOST_ROOT) + endif + CXXFLAGS += -DWARPX_QED_TABLE_GEN + CFLAGS += -DWARPX_QED_TABLE_GEN + FFLAGS += -DWARPX_QED_TABLE_GEN + F90FLAGS += -DWARPX_QED_TABLE_GEN + USERSuffix := $(USERSuffix).GENTABLES + endif +endif include $(PICSAR_HOME)/src/Make.package diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 3341a8fe8..f9a0e51d7 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -19,6 +19,8 @@ #include <map> #include <string> #include <algorithm> +#include <utility> +#include <tuple> // // MultiParticleContainer holds multiple (nspecies or npsecies+1 when @@ -218,13 +220,58 @@ protected: #ifdef WARPX_QED // The QED engines - std::shared_ptr<BreitWheelerEngine> shr_p_bw_engine; - std::shared_ptr<QuantumSynchrotronEngine> shr_p_qs_engine; + std::shared_ptr<BreitWheelerEngine> m_shr_p_bw_engine; + std::shared_ptr<QuantumSynchrotronEngine> m_shr_p_qs_engine; //_______________________________ - //Initialize QED engines and provides smart pointers - //to species who need QED processes + /** + * Initialize QED engines and provides smart pointers + * to species who need QED processes + */ void InitQED (); + + //Variables to store how many species need a QED process + int m_nspecies_quantum_sync = 0; + int m_nspecies_breit_wheeler = 0; + //________ + + /** + * Returns the number of species having Quantum Synchrotron process enabled + */ + int NSpeciesQuantumSync() const { return m_nspecies_quantum_sync;} + + /** + * Returns the number of species having Breit Wheeler process enabled + */ + int NSpeciesBreitWheeler() const { return m_nspecies_breit_wheeler;} + + /** + * Initializes the Quantum Synchrotron engine + */ + void InitQuantumSync (); + + /** + * Initializes the Quantum Synchrotron engine + */ + void InitBreitWheeler (); + + /** + * Parses inputfile parameters for Quantum Synchrotron engine + * @return {a tuple containing a flag which is true if tables + * have to be generate, a filename (where tables should be stored + * or read from) and control parameters.} + */ + std::tuple<bool, std::string, PicsarQuantumSynchrotronCtrl> + ParseQuantumSyncParams (); + + /** + * Parses inputfile parameters for Breit Wheeler engine + * @return {a tuple containing a flag which is true if tables + * have to be generate, a filename (where tables should be stored + * or read from) and control parameters.} + */ + std::tuple<bool, std::string, PicsarBreitWheelerCtrl> + ParseBreitWheelerParams (); #endif private: diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index de2b4583d..fca22daa2 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,11 +1,16 @@ #include <MultiParticleContainer.H> + +#include <AMReX_Vector.H> + #include <WarpX_f.H> #include <WarpX.H> +//This is now needed for writing a binary file on disk. +#include <WarpXUtil.H> + #include <limits> #include <algorithm> #include <string> -#include <memory> using namespace amrex; @@ -615,18 +620,295 @@ MultiParticleContainer::doFieldIonization () #ifdef WARPX_QED void MultiParticleContainer::InitQED () { - shr_p_qs_engine = std::make_shared<QuantumSynchrotronEngine>(); - shr_p_bw_engine = std::make_shared<BreitWheelerEngine>(); + m_shr_p_qs_engine = std::make_shared<QuantumSynchrotronEngine>(); + m_shr_p_bw_engine = std::make_shared<BreitWheelerEngine>(); + + m_nspecies_quantum_sync = 0; + m_nspecies_breit_wheeler = 0; for (auto& pc : allcontainers) { if(pc->has_quantum_sync()){ pc->set_quantum_sync_engine_ptr - (shr_p_qs_engine); + (m_shr_p_qs_engine); + m_nspecies_quantum_sync++; } if(pc->has_breit_wheeler()){ pc->set_breit_wheeler_engine_ptr - (shr_p_bw_engine); + (m_shr_p_bw_engine); + m_nspecies_breit_wheeler++; } } + + if(m_nspecies_quantum_sync != 0) + InitQuantumSync(); + + if(m_nspecies_breit_wheeler !=0) + InitBreitWheeler(); + +} + +void MultiParticleContainer::InitQuantumSync () +{ + bool generate_table; + PicsarQuantumSynchrotronCtrl ctrl; + std::string filename; + std::tie(generate_table, filename, ctrl) = ParseQuantumSyncParams(); + + //Only temporary for test purposes, will be removed + ParmParse pp("qed_qs"); + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) return; + //_________________________________________________ + + + if(generate_table && ParallelDescriptor::IOProcessor()){ + m_shr_p_qs_engine->compute_lookup_tables(ctrl); + Vector<char> all_data = m_shr_p_qs_engine->export_lookup_tables_data(); + WarpXUtilIO::WriteBinaryDataOnFile(filename, all_data); + } + ParallelDescriptor::Barrier(); + + Vector<char> table_data; + ParallelDescriptor::ReadAndBcastFile(filename, table_data); + ParallelDescriptor::Barrier(); + + //No need to initialize from raw data for the processor that + //has just generated the table + if(!generate_table || !ParallelDescriptor::IOProcessor()){ + m_shr_p_qs_engine->init_lookup_tables_from_raw_data(table_data); + } + + if(!m_shr_p_qs_engine->are_lookup_tables_initialized()) + amrex::Error("Table initialization has failed!\n"); +} + +void MultiParticleContainer::InitBreitWheeler () +{ + bool generate_table; + PicsarBreitWheelerCtrl ctrl; + std::string filename; + std::tie(generate_table, filename, ctrl) = ParseBreitWheelerParams(); + + //Only temporary for test purposes, will be removed + ParmParse pp("qed_bw"); + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) return; + //_________________________________________________ + + if(generate_table && ParallelDescriptor::IOProcessor()){ + m_shr_p_bw_engine->compute_lookup_tables(ctrl); + Vector<char> all_data = m_shr_p_bw_engine->export_lookup_tables_data(); + WarpXUtilIO::WriteBinaryDataOnFile(filename, all_data); + } + ParallelDescriptor::Barrier(); + + Vector<char> table_data; + ParallelDescriptor::ReadAndBcastFile(filename, table_data); + ParallelDescriptor::Barrier(); + + //No need to initialize from raw data for the processor that + //has just generated the table + if(!generate_table || !ParallelDescriptor::IOProcessor()){ + m_shr_p_bw_engine->init_lookup_tables_from_raw_data(table_data); + } + + if(!m_shr_p_bw_engine->are_lookup_tables_initialized()) + amrex::Error("Table initialization has failed!\n"); +} + +std::tuple<bool,std::string,PicsarQuantumSynchrotronCtrl> +MultiParticleContainer::ParseQuantumSyncParams () +{ + PicsarQuantumSynchrotronCtrl ctrl = + m_shr_p_qs_engine->get_default_ctrl(); + bool generate_table{false}; + std::string table_name; + + ParmParse pp("qed_qs"); + + // Engine paramenter: chi_part_min is the minium chi parameter to be + // considered by the engine. If a lepton has chi < chi_part_min, + // the optical depth is not evolved and photon generation is ignored + pp.query("chi_min", ctrl.chi_part_min); + + //Only temporary for test purposes, will be removed + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) + return std::make_tuple(false, "__DUMMY__", ctrl); + //_________________________________________________ + + pp.query("generate_table", generate_table); + if(generate_table){ + int t_int = 0; + + //==Table parameters== + + //--- sub-table 1 (1D) + //These parameters are used to pre-compute a function + //which appears in the evolution of the optical depth + + //Minimun chi for the table. If a lepton has chi < chi_part_tdndt_min, + ///chi is considered as it were equal to chi_part_tdndt_min + pp.query("tab_dndt_chi_min", ctrl.chi_part_tdndt_min); + + //Maximum chi for the table. If a lepton has chi > chi_part_tdndt_max, + ///chi is considered as it were equal to chi_part_tdndt_max + pp.query("tab_dndt_chi_max", ctrl.chi_part_tdndt_max); + + //How many points should be used for chi in the table + pp.query("tab_dndt_how_many", t_int); + ctrl.chi_part_tdndt_how_many= t_int; + //------ + + //--- sub-table 2 (2D) + //These parameters are used to pre-compute a function + //which is used to extract the properties of the generated + //photons. + + //Minimun chi for the table. If a lepton has chi < chi_part_tem_min, + ///chi is considered as it were equal to chi_part_tem_min + pp.query("tab_em_chi_min", ctrl.chi_part_tem_min); + + //Maximum chi for the table. If a lepton has chi > chi_part_tem_max, + ///chi is considered as it were equal to chi_part_tem_max + pp.query("tab_em_chi_max", ctrl.chi_part_tem_max); + + //How many points should be used for chi in the table + pp.query("tab_em_chi_how_many", t_int); + ctrl.chi_part_tem_how_many = t_int; + + //The other axis of the table is a cumulative probability distribution + //(corresponding to different energies of the generated particles) + //This parameter is the number of different points to consider + pp.query("tab_em_prob_how_many", t_int); + ctrl.prob_tem_how_many = t_int; + //==================== + + pp.query("save_table_in", table_name); + + } + + std::string load_table_name; + pp.query("load_table_from", load_table_name); + if(!load_table_name.empty()){ + if(generate_table && ParallelDescriptor::IOProcessor()){ + amrex::Print() << "Warning, Quantum Synchrotron table will be loaded, not generated. \n"; + } + table_name = load_table_name; + generate_table = false; + } + +#ifndef WARPX_QED_TABLE_GEN + if(generate_table){ + amrex::Error("Error: Compile with QED_TABLE_GEN=TRUE to enable table generation!\n"); + } +#endif + + if(table_name.empty()){ + amrex::Error("Error: Quantum Synchrotron table has either to be generated or to be loaded.\n"); + } + + return std::make_tuple(generate_table, table_name, ctrl); +} + +std::tuple<bool,std::string,PicsarBreitWheelerCtrl> +MultiParticleContainer::ParseBreitWheelerParams () +{ + PicsarBreitWheelerCtrl ctrl = + m_shr_p_bw_engine->get_default_ctrl(); + bool generate_table{false}; + std::string table_name; + + ParmParse pp("qed_bw"); + + // Engine paramenter: chi_phot_min is the minium chi parameter to be + // considered by the engine. If a photon has chi < chi_phot_min, + // the optical depth is not evolved and pair generation is ignored + pp.query("chi_min", ctrl.chi_phot_min); + + //Only temporary for test purposes, will be removed + bool ignore_tables = false; + pp.query("ignore_tables_for_test", ignore_tables); + if(ignore_tables) + return std::make_tuple(false, "__DUMMY__", ctrl); + //_________________________________________________ + + pp.query("generate_table", generate_table); + if(generate_table){ + + int t_int; + + //==Table parameters== + + //--- sub-table 1 (1D) + //These parameters are used to pre-compute a function + //which appears in the evolution of the optical depth + + //Minimun chi for the table. If a photon has chi < chi_phot_tdndt_min, + //an analytical approximation is used. + pp.query("tab_dndt_chi_min", ctrl.chi_phot_tdndt_min); + + //Maximum chi for the table. If a photon has chi > chi_phot_tdndt_min, + //an analytical approximation is used. + pp.query("tab_dndt_chi_max", ctrl.chi_phot_tdndt_max); + + //How many points should be used for chi in the table + pp.query("tab_dndt_how_many", t_int); + ctrl.chi_phot_tdndt_how_many = t_int; + //------ + + //--- sub-table 2 (2D) + //These parameters are used to pre-compute a function + //which is used to extract the properties of the generated + //particles. + + //Minimun chi for the table. If a photon has chi < chi_phot_tpair_min + //chi is considered as it were equal to chi_phot_tpair_min + pp.query("tab_pair_chi_min", ctrl.chi_phot_tpair_min); + + //Maximum chi for the table. If a photon has chi > chi_phot_tpair_max + //chi is considered as it were equal to chi_phot_tpair_max + pp.query("tab_pair_chi_max", ctrl.chi_phot_tpair_max); + + //How many points should be used for chi in the table + pp.query("tab_pair_chi_how_many", t_int); + ctrl.chi_phot_tpair_how_many = t_int; + + //The other axis of the table is the fraction of the initial energy + //'taken away' by the most energetic particle of the pair. + //This parameter is the number of different fractions to consider + pp.query("tab_pair_frac_how_many", t_int); + ctrl.chi_frac_tpair_how_many = t_int; + //==================== + + pp.query("save_table_in", table_name); + } + + std::string load_table_name; + pp.query("load_table_from", load_table_name); + if(!load_table_name.empty()){ + if(generate_table && ParallelDescriptor::IOProcessor()){ + amrex::Print() << "Warning, Breit Wheeler table will be loaded, not generated. \n"; + } + table_name = load_table_name; + generate_table = false; + } + +#ifndef WARPX_QED_TABLE_GEN + if(generate_table){ + if(ParallelDescriptor::IOProcessor()){ + amrex::Error("Error: Compile with QED_TABLE_GEN=TRUE to enable table generation!\n"); + } + } +#endif + + if(table_name.empty()){ + amrex::Error("Error: Breit Wheeler table has either to be generated or to be loaded.\n"); + } + + return std::make_tuple(generate_table, table_name, ctrl); } #endif diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 612da01ca..7e52b52e1 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -26,9 +26,9 @@ PhotonParticleContainer::PhotonParticleContainer (AmrCore* amr_core, int ispecie ParmParse pp(species_name); #ifdef WARPX_QED - //IF do_qed is enabled, find out if Breit Wheeler process is enabled - if(do_qed) - pp.query("do_qed_breit_wheeler", do_qed_breit_wheeler); + //IF m_do_qed is enabled, find out if Breit Wheeler process is enabled + if(m_do_qed) + pp.query("do_qed_breit_wheeler", m_do_qed_breit_wheeler); //Check for processes which do not make sense for photons bool test_quantum_sync = false; diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b18c9b5f8..17a504719 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> @@ -219,16 +220,16 @@ protected: #ifdef WARPX_QED // A flag to enable quantum_synchrotron process for leptons - bool do_qed_quantum_sync = false; + bool m_do_qed_quantum_sync = false; // A flag to enable breit_wheeler process [photons only!!] - bool do_qed_breit_wheeler = false; + bool m_do_qed_breit_wheeler = false; // A smart pointer to an instance of a Quantum Synchrotron engine - std::shared_ptr<QuantumSynchrotronEngine> shr_ptr_qs_engine; + std::shared_ptr<QuantumSynchrotronEngine> m_shr_p_qs_engine; // A smart pointer to an instance of a Breit Wheeler engine [photons only!] - std::shared_ptr<BreitWheelerEngine> shr_ptr_bw_engine; + std::shared_ptr<BreitWheelerEngine> m_shr_p_bw_engine; #endif }; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index fed7266e1..51690d659 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -74,13 +74,13 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp #ifdef WARPX_QED //Add real component if QED is enabled - pp.query("do_qed", do_qed); - if(do_qed) + pp.query("do_qed", m_do_qed); + if(m_do_qed) AddRealComp("tau"); //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled - if(do_qed) - pp.query("do_qed_quantum_sync", do_qed_quantum_sync); + if(m_do_qed) + pp.query("do_qed_quantum_sync", m_do_qed_quantum_sync); //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!! #endif @@ -91,7 +91,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp plot_flag_size += 6; #ifdef WARPX_QED - if(do_qed){ + if(m_do_qed){ // plot_flag will have an entry for the optical depth plot_flag_size++; } @@ -123,7 +123,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } #ifdef WARPX_QED - if(do_qed){ + if(m_do_qed){ //Optical depths is always plotted if QED is on plot_flags[plot_flag_size-1] = 1; } @@ -544,11 +544,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) BreitWheelerGetOpticalDepth breit_wheeler_get_opt; if(loc_has_quantum_sync){ quantum_sync_get_opt = - shr_ptr_qs_engine->build_optical_depth_functor(); + m_shr_p_qs_engine->build_optical_depth_functor(); } if(loc_has_breit_wheeler){ breit_wheeler_get_opt = - shr_ptr_bw_engine->build_optical_depth_functor(); + m_shr_p_bw_engine->build_optical_depth_functor(); } #endif @@ -2238,25 +2238,25 @@ PhysicalParticleContainer::AmIALepton(){ bool PhysicalParticleContainer::has_quantum_sync() { - return do_qed_quantum_sync; + return m_do_qed_quantum_sync; } bool PhysicalParticleContainer::has_breit_wheeler() { - return do_qed_breit_wheeler; + return m_do_qed_breit_wheeler; } void PhysicalParticleContainer:: set_breit_wheeler_engine_ptr(std::shared_ptr<BreitWheelerEngine> ptr) { - shr_ptr_bw_engine = ptr; + m_shr_p_bw_engine = ptr; } void PhysicalParticleContainer:: set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine> ptr) { - shr_ptr_qs_engine = ptr; + m_shr_p_qs_engine = ptr; } #endif diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index b41dcede3..dbd913c5b 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -320,7 +320,7 @@ protected: int do_back_transformed_diagnostics = 1; #ifdef WARPX_QED - bool do_qed = false; + bool m_do_qed = false; virtual bool has_quantum_sync(){return false;}; virtual bool has_breit_wheeler(){return false;}; diff --git a/Source/QED/BreitWheelerEngineInnards.H b/Source/QED/BreitWheelerEngineInnards.H new file mode 100644 index 000000000..640cdfa94 --- /dev/null +++ b/Source/QED/BreitWheelerEngineInnards.H @@ -0,0 +1,46 @@ +#ifndef WARPX_breit_wheeler_engine_innards_h_ +#define WARPX_breit_wheeler_engine_innards_h_ + +#include "QedWrapperCommons.H" + +#include <AMReX_Gpu.H> + +//This includes only the definition of a simple datastructure +//used to control the Breit Wheeler engine. +#include <breit_wheeler_engine_ctrl.h> + +/** + * This structure holds all the parameters required to use the + * Breit Wheeler engine: a POD control structure and lookup + * tables data. + */ +struct BreitWheelerEngineInnards +{ + // Control parameters (a POD struct) + // ctrl contains several parameters: + // - chi_phot_min : the minium chi parameter to be + // considered by the engine + // - chi_phot_tdndt_min : minimun chi for sub-table 1 (1D) + // - chi_phot_tdndt_max : maximum chi for sub-table 1 (1D) + // - chi_phot_tdndt_how_many : how many points to use for sub-table 1 (1D) + // - chi_phot_tpair_min : minimun chi for sub-table 2 (1D) + // - chi_phot_tpair_max : maximum chi for sub-table 2 (1D) + // - chi_phot_tpair_how_many : how many points to use for chi for sub-table 2 (2D) + // - chi_frac_tpair_how_many : how many points to use for the second axis of sub-table 2 (2D) + picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real> ctrl; + + //Lookup table data + //---sub-table 1 (1D) + amrex::Gpu::ManagedVector<amrex::Real> TTfunc_coords; + amrex::Gpu::ManagedVector<amrex::Real> TTfunc_data; + //--- + + //---sub-table 2 (2D) + 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; + //______ +}; +//========================================================== + +#endif //WARPX_breit_wheeler_engine_innards_h_ diff --git a/Source/QED/BreitWheelerEngineTableBuilder.H b/Source/QED/BreitWheelerEngineTableBuilder.H new file mode 100644 index 000000000..e30b82208 --- /dev/null +++ b/Source/QED/BreitWheelerEngineTableBuilder.H @@ -0,0 +1,26 @@ +#ifndef WARPX_breit_wheeler_engine_table_builder_h_ +#define WARPX_breit_wheeler_engine_table_builder_h_ + +#include "QedWrapperCommons.H" +#include "BreitWheelerEngineInnards.H" + +//This includes only the definition of a simple datastructure +//used to control the Breit Wheeler engine. +#include <breit_wheeler_engine_ctrl.h> + +/** + * A class which computes the lookup tables for the Breit Wheeler engine. + */ +class BreitWheelerEngineTableBuilder{ + public: + /** + * Computes the tables. + * @param[in] ctrl control parameters to generate the tables + * @param[out] innards structure holding both a copy of ctrl and lookup tables data + */ + void compute_table + (picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real> ctrl, + BreitWheelerEngineInnards& innards) const; +}; + +#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..3326d5b59 --- /dev/null +++ b/Source/QED/BreitWheelerEngineTableBuilder.cpp @@ -0,0 +1,58 @@ +#include "BreitWheelerEngineTableBuilder.H" + +//Include the full Breit Wheeler engine with table generation support +//(after some consistency tests). This requires to have a recent version +// of the Boost library. +#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> +//_______________________________________________ + +//Some handy aliases +using PicsarBreitWheelerEngine = picsar::multi_physics:: + breit_wheeler_engine<amrex::Real, QedUtils::DummyStruct>; + +using PicsarBreitWheelerCtrl = + picsar::multi_physics::breit_wheeler_engine_ctrl<amrex::Real>; +//_______________________________________________ + +void +BreitWheelerEngineTableBuilder::compute_table + (PicsarBreitWheelerCtrl ctrl, + BreitWheelerEngineInnards& innards) const +{ + 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); + innards.cum_distrib_coords_1.assign( + bw_innards_picsar.cum_distrib_table_coords_1_ptr, + bw_innards_picsar.cum_distrib_table_coords_1_ptr + + bw_innards_picsar.cum_distrib_table_coords_1_how_many); + innards.cum_distrib_coords_2.assign( + bw_innards_picsar.cum_distrib_table_coords_2_ptr, + bw_innards_picsar.cum_distrib_table_coords_2_ptr + + bw_innards_picsar.cum_distrib_table_coords_2_how_many); + innards.cum_distrib_data.assign( + bw_innards_picsar.cum_distrib_table_data_ptr, + bw_innards_picsar.cum_distrib_table_data_ptr + + bw_innards_picsar.cum_distrib_table_data_how_many); + //____ +} diff --git a/Source/QED/BreitWheelerEngineWrapper.H b/Source/QED/BreitWheelerEngineWrapper.H index 7cdc66b44..1033ff7c9 100644 --- a/Source/QED/BreitWheelerEngineWrapper.H +++ b/Source/QED/BreitWheelerEngineWrapper.H @@ -2,53 +2,302 @@ #define WARPX_breit_wheeler_engine_wrapper_h_ #include "QedWrapperCommons.H" +#include "BreitWheelerEngineInnards.H" -//BW ENGINE from PICSAR +#include <AMReX_Array.H> +#include <AMReX_Vector.H> +#include <AMReX_Gpu.H> + +//#define PXRMP_CORE_ONLY allows importing only the 'core functions' of the +//Breit Wheeler engine of the QED PICSAR library. #define PXRMP_CORE_ONLY -#include "breit_wheeler_engine.hpp" +#include <breit_wheeler_engine.hpp> + +//Lookup table building function is in a dedicated (optional) class to +//avoid including heavy dependencies if they are not needed. +#ifdef WARPX_QED_TABLE_GEN + #include "BreitWheelerEngineTableBuilder.H" +#endif + +#include <string> + +//Some handy aliases + +// The engine has two templated arguments: the numerical type +// and a random number generator. However, random numbers are not +// used to generate the lookup tables and the static member +// functions which are called in the functors do not use +// random numbers as well. Therefore, an empty "DummyStruct" +// can be passed as a template parameter. +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>; +//__________ // Functors ================================== -// These functors provide the core elementary functions of the library -// Can be included in GPU kernels +// These functors allow using the core elementary functions of the library. +// They are generated by a factory class (BreitWheelerEngine, see below). +// They can be included in GPU kernels. /** - * \brief Functor to initialize the optical depth of photons for the + * Functor to initialize the optical depth of photons for the * Breit-Wheeler process */ class BreitWheelerGetOpticalDepth { public: + /** + * Constructor does nothing because optical depth initialization + * does not require control parameters or lookup tables. + */ BreitWheelerGetOpticalDepth () {}; + /** + * () operator is just a thin wrapper around a very simple function to + * generate the optical depth. It can be used on GPU. + */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE - amrex::Real operator() () const + amrex::Real operator() () const noexcept { - return WarpXBreitWheelerWrapper:: + //A random number in [0,1) should be provided as an argument. + return PicsarBreitWheelerEngine:: internal_get_optical_depth(amrex::Random()); } }; //____________________________________________ +/** + * Functor to evolve the optical depth of photons due to the + * Breit-Wheeler process + */ +class BreitWheelerEvolveOpticalDepth +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + BreitWheelerEvolveOpticalDepth(BreitWheelerEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_TTfunc_size{r_innards.TTfunc_coords.size()}, + m_p_TTfunc_coords{r_innards.TTfunc_coords.dataPtr()}, + m_p_TTfunc_data{r_innards.TTfunc_data.dataPtr()} + {}; + + /** + * Evolves the optical depth. It can be used on GPU. + * @param[in] px,py,pz momentum components of the photon (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] dt timestep (SI units) + * @param[in,out] opt_depth optical depth of the photon. It is modified by the method. + * @return a flag which is 1 if optical depth becomes negative (i.e. a pair has to be generated). + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + int 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 noexcept + { + bool has_event_happened{false}; + + //the library provides the time (< dt) at which the event occurs, but this + //feature won't be used in WarpX for now. + amrex::Real unused_event_time{0.0}; + + PicsarBreitWheelerEngine:: + internal_evolve_opt_depth_and_determine_event( + px, py, pz, + ex, ey, ez, + bx, by, bz, + dt, opt_depth, + has_event_happened, unused_event_time, + m_dummy_lambda, + picsar::multi_physics::lookup_1d<amrex::Real>{ + m_TTfunc_size, + m_p_TTfunc_coords, + m_p_TTfunc_data}, + m_ctrl); + + return has_event_happened; + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarBreitWheelerCtrl m_ctrl; + + //lookup table data + size_t m_TTfunc_size; + amrex::Real* m_p_TTfunc_coords; + amrex::Real* m_p_TTfunc_data; +}; + +/** + * Functor to generate a pair via the + * Breit-Wheeler process + */ +class BreitWheelerGeneratePairs +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + BreitWheelerGeneratePairs( + BreitWheelerEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_cum_distrib_coords_1_size{r_innards.cum_distrib_coords_1.size()}, + m_cum_distrib_coords_2_size{r_innards.cum_distrib_coords_2.size()}, + m_p_distrib_coords_1{r_innards.cum_distrib_coords_1.data()}, + m_p_distrib_coords_2{r_innards.cum_distrib_coords_2.data()}, + m_p_cum_distrib_data{r_innards.cum_distrib_data.data() + }{}; + + /** + * Generates sampling (template parameter) pairs according to Breit Wheeler process. + * It can be used on GPU. + * @param[in] px,py,pz momentum components of the photon (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] weight of the photon (code units) + * @param[out] e_px,e_py,e_pz momenta of generated electrons. Each array should have size=sampling (SI units) + * @param[out] p_px,p_py,p_pz momenta of generated positrons. Each array should have size=sampling (SI units) + * @param[out] e_weight,p_weight weight of the generated particles Each array should have size=sampling (code units). + */ + template <size_t sampling> + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void 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 noexcept + { + //[sampling] random numbers are needed + picsar::multi_physics::picsar_array<amrex::Real, sampling> + rand_zero_one_minus_epsi; + for(auto& el : rand_zero_one_minus_epsi) el = amrex::Random(); + + 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, + m_dummy_lambda, + picsar::multi_physics::lookup_2d<amrex::Real>{ + m_cum_distrib_coords_1_size, + m_p_distrib_coords_1, + m_cum_distrib_coords_2_size, + m_p_distrib_coords_2, + m_p_cum_distrib_data + }, + m_ctrl, + rand_zero_one_minus_epsi.data()); + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarBreitWheelerCtrl m_ctrl; + + //lookup table data + size_t m_cum_distrib_coords_1_size; + size_t m_cum_distrib_coords_2_size; + amrex::Real* m_p_distrib_coords_1; + amrex::Real* m_p_distrib_coords_2; + amrex::Real* m_p_cum_distrib_data; +}; + // Factory class ============================= /** - * \brief Wrapper for the Breit Wheeler engine of the PICSAR library + * Wrapper for the Breit Wheeler engine of the PICSAR library */ class BreitWheelerEngine { public: + /** + * Constructor requires no arguments. + */ BreitWheelerEngine (); /** - * \brief Builds the functor to initialize the optical depth + * Builds the functor to initialize the optical depth */ BreitWheelerGetOpticalDepth build_optical_depth_functor (); + + /** + * Builds the functor to evolve the optical depth + */ + BreitWheelerEvolveOpticalDepth build_evolve_functor (); + + /** + * Builds the functor to generate the pairs + */ + BreitWheelerGeneratePairs build_pair_functor (); + + /** + * Checks if the optical tables are properly initialized + */ + bool are_lookup_tables_initialized () const; + + /** + * Init lookup tables from raw binary data. + * @param[in] raw_data a Vector of char + * @return true if it succeeds, false if it cannot parse raw_data + */ + bool init_lookup_tables_from_raw_data (const amrex::Vector<char>& raw_data); + + /** + * Export lookup tables data into a raw binary Vector + * @return the data in binary format. The Vector is empty if tables were + * not previously initialized. + */ + amrex::Vector<char> export_lookup_tables_data () const; + + /** + * Computes the lookup tables. It does nothing unless WarpX is compiled with QED_TABLE_GEN=TRUE + * @param[in] ctrl control params to generate the tables + */ + void compute_lookup_tables (PicsarBreitWheelerCtrl ctrl); + + /** + * gets default (reasonable) values for the control parameters + * @return default control params to generate the tables + */ + PicsarBreitWheelerCtrl get_default_ctrl() const; + +private: + bool m_lookup_tables_initialized = false; + + BreitWheelerEngineInnards m_innards; + +//Table builing is available only if WarpX is compiled with QED_TABLE_GEN=TRUE +#ifdef WARPX_QED_TABLE_GEN + BreitWheelerEngineTableBuilder m_table_builder; +#endif + }; //============================================ diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp index 33e6d48d7..42953c97f 100644 --- a/Source/QED/BreitWheelerEngineWrapper.cpp +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -1,16 +1,188 @@ #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 () +BreitWheelerGetOpticalDepth +BreitWheelerEngine::build_optical_depth_functor () { return BreitWheelerGetOpticalDepth(); } +BreitWheelerEvolveOpticalDepth +BreitWheelerEngine::build_evolve_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return BreitWheelerEvolveOpticalDepth(m_innards); +} + +BreitWheelerGeneratePairs +BreitWheelerEngine::build_pair_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return BreitWheelerGeneratePairs(m_innards); +} + +bool BreitWheelerEngine::are_lookup_tables_initialized () const +{ + return m_lookup_tables_initialized; +} + +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, m_innards.ctrl.chi_phot_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tdndt_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tdndt_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tdndt_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tdndt_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tdndt_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tdndt_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tpair_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tpair_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tpair_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tpair_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_phot_tpair_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_phot_tpair_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_frac_tpair_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_frac_tpair_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + //___________________________ + + //Data + Vector<Real> tndt_coords(m_innards.ctrl.chi_phot_tdndt_how_many); + Vector<Real> tndt_data(m_innards.ctrl.chi_phot_tdndt_how_many); + Vector<Real> cum_tab_coords1(m_innards.ctrl.chi_phot_tpair_how_many); + Vector<Real> cum_tab_coords2(m_innards.ctrl.chi_frac_tpair_how_many); + Vector<Real> cum_tab_data(m_innards.ctrl.chi_phot_tpair_how_many* + m_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; + m_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; + m_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; + m_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; + m_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; + m_innards.cum_distrib_data.assign( + cum_tab_data.begin(), cum_tab_data.end()); + + //___________________________ + m_lookup_tables_initialized = true; + + return true; +} + +Vector<char> BreitWheelerEngine::export_lookup_tables_data () const +{ + Vector<char> res{}; + + if(!m_lookup_tables_initialized) + return res; + + add_data_to_vector_char(&m_innards.ctrl.chi_phot_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tdndt_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tdndt_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tdndt_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tpair_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tpair_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_phot_tpair_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_frac_tpair_how_many, 1, res); + + add_data_to_vector_char(m_innards.TTfunc_coords.data(), + m_innards.TTfunc_coords.size(), res); + add_data_to_vector_char(m_innards.TTfunc_data.data(), + m_innards.TTfunc_data.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_1.data(), + m_innards.cum_distrib_coords_1.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_2.data(), + m_innards.cum_distrib_coords_2.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_data.data(), + m_innards.cum_distrib_data.size(), res); + + return res; +} + +PicsarBreitWheelerCtrl +BreitWheelerEngine::get_default_ctrl() const +{ + return PicsarBreitWheelerCtrl(); +} + +void BreitWheelerEngine::compute_lookup_tables ( + PicsarBreitWheelerCtrl ctrl) +{ +#ifdef WARPX_QED_TABLE_GEN + m_table_builder.compute_table(ctrl, m_innards); + m_lookup_tables_initialized = true; +#endif +} + //============================================ diff --git a/Source/QED/Make.package b/Source/QED/Make.package index d4bad3f18..c9cd73cc2 100644 --- a/Source/QED/Make.package +++ b/Source/QED/Make.package @@ -1,8 +1,21 @@ CEXE_headers += QedWrapperCommons.H +CEXE_headers += QedChiFunctions.H +CEXE_headers += QedTableParserHelperFunctions.H +CEXE_headers += BreitWheelerEngineInnards.H +CEXE_headers += QuantumSyncEngineInnards.H CEXE_headers += BreitWheelerEngineWrapper.H -CEXE_headers += QuantumSyncsEngineWrapper.H +CEXE_headers += QuantumSyncEngineWrapper.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..dd8ffac0e --- /dev/null +++ b/Source/QED/QedChiFunctions.H @@ -0,0 +1,62 @@ +#ifndef WARPX_amrex_qed_chi_functions_h_ +#define WARPX_amrex_qed_chi_functions_h_ + +/** + * This header contains wrappers around functions provided by + * the PICSAR QED library to calculate the 'chi' parameter + * for photons or electrons and positrons. + */ + +#include "QedWrapperCommons.H" + +//#define PXRMP_CORE_ONLY allows importing only the 'core functions' of the +//QED library. +#define PXRMP_CORE_ONLY +#include <chi_functions.hpp> + +namespace QedUtils{ + /** + * Function to calculate the 'chi' parameter for photons. + * Suitable for GPU kernels. + * @param[in] px,py,pz components of momentum (SI units) + * @param[in] ex,ey,ez components of electric field (SI units) + * @param[in] bx,by,bz components of magnetic field (SI units) + * @return chi parameter + */ + 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) + { + //laser wavelength is unused if SI units are set + const amrex::Real dummy_lambda = 1.0; + return picsar::multi_physics::chi_photon( + px, py, pz, ex, ey, ez, bx, by, bz, dummy_lambda); + } + + /** + * Function to calculate the 'chi' parameter for electrons or positrons. + * Suitable for GPU kernels. + * @param[in] px,py,pz components of momentum (SI units) + * @param[in] ex,ey,ez components of electric field (SI units) + * @param[in] bx,by,bz components of magnetic field (SI units) + * @return chi parameter + */ + 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) + { + //laser wavelength is unused if SI units are set + 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_ diff --git a/Source/QED/QedTableParserHelperFunctions.H b/Source/QED/QedTableParserHelperFunctions.H new file mode 100644 index 000000000..9f9f37017 --- /dev/null +++ b/Source/QED/QedTableParserHelperFunctions.H @@ -0,0 +1,90 @@ +#ifndef WARPX_amrex_qed_table_parser_helper_functions_h_ +#define WARPX_amrex_qed_table_parser_helper_functions_h_ + +/** + * This header contains helper functions to safely extract data + * (e.g. integers, floating point numbers) from raw binary data + * (i.e. a char*) and to convert arrays into raw binary data. + */ + +#include <AMReX_Vector.H> +#include <tuple> + +namespace QedUtils{ + /** + * This function safely extracts an amrex::Vector<T> from raw binary data. + * T must be a simple datatype (e.g. an int, a float, a double...). + * + * @param[in] p_char a pointer to the binary stream + * @param[in] how_many how many T should be read from stream + * @param[in] p_last a pointer to the last element of the char* array + * @return {a tuple containing + * 1) flag (which is false if p_last is exceeded) + * 2) a Vector of T + * 3) a pointer to a new location of the binary data (after having read how_many T)} + */ + template <class T> + std::tuple<bool, amrex::Vector<T>, const char*>parse_raw_data_vec( + const char* p_data, size_t how_many, const char* const p_last) + { + amrex::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); + } + + /** + * This function safely extracts a T from raw binary data. + * T must be a simple datatype (e.g. an int, a float, a double...). + * + * @param[in] p_char a pointer to the binary stream + * @param[in] p_last a pointer to the last element of the char* array + * @return {a tuple containing + * 1) flag (which is false if p_last is exceeded) + * 2) a T + * 3) a pointer to a new location of the binary data (after having read 1 T)} + */ + 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); + } + + /** + * This function converts a C-style array of T into + * a Vector<char> (i.e. raw binary data) and adds it + * to an existing Vector<char> passed by reference + * @param[in] p_data a pointer to the beginning of the array + * @param[in] how_many number of elements of type T in the array + * @param[in,out] raw_data data will be appended to this vector + */ + template <class T> + void add_data_to_vector_char ( + const T* p_data, size_t how_many, amrex::Vector<char>& raw_data) + { + raw_data.insert( + raw_data.end(), + reinterpret_cast<const char*>(p_data), + reinterpret_cast<const char*>(p_data) + + sizeof(T)*how_many + ); + } +}; + +#endif //WARPX_amrex_qed_table_parser_helper_functions_h_ diff --git a/Source/QED/QedWrapperCommons.H b/Source/QED/QedWrapperCommons.H index 821034c06..210e7247e 100644 --- a/Source/QED/QedWrapperCommons.H +++ b/Source/QED/QedWrapperCommons.H @@ -1,16 +1,36 @@ #ifndef WARPX_amrex_qed_wrapper_commons_h_ #define WARPX_amrex_qed_wrapper_commons_h_ -//Common definitions for the QED library wrappers +/** + * This header contains some common #define directives and a + * 'dummy' class used by the QED library wrappers and related + * components. + */ #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 +/** + * PICSAR uses PXRMP_GPU to decorate methods which should be + * compiled for GPU. The user has to set it to the right value + * (AMREX_GPU_DEVICE in this case). + * PXRMP_WITH_SI_UNITS sets the library to use International + * System units. + */ +#define PXRMP_GPU AMREX_GPU_HOST_DEVICE #define PXRMP_WITH_SI_UNITS +//_________________________ + +/** + * A namespace called 'QedUtils' is used to encapsulate + * free functions (defined elsewhere) and an + * empty datastructure (DummyStruct), which is re-used by several + * components. + */ +namespace QedUtils{ + struct DummyStruct{}; +}; +//_________________________ -//An empty data type -struct DummyStruct{}; #endif //WARPX_amrex_qed_wrapper_commons_h_ diff --git a/Source/QED/QuantumSyncEngineInnards.H b/Source/QED/QuantumSyncEngineInnards.H new file mode 100644 index 000000000..6206b103a --- /dev/null +++ b/Source/QED/QuantumSyncEngineInnards.H @@ -0,0 +1,46 @@ +#ifndef WARPX_quantum_sync_engine_innards_h_ +#define WARPX_quantum_sync_engine_innards_h_ + +#include "QedWrapperCommons.H" + +#include <AMReX_Gpu.H> + +//This includes only the definition of a simple datastructure +//used to control the Quantum Synchrotron engine. +#include <quantum_sync_engine_ctrl.h> + +/** + * This structure holds all the parameters required to use the + * Quantum Synchrotron engine: a POD control structure and lookup + * tables data. + */ +struct QuantumSynchrotronEngineInnards +{ + // Control parameters (a POD struct) + // ctrl contains several parameters: + // - chi_part_min : the minium chi parameter to be + // considered by the engine + // - chi_part_tdndt_min : minimun chi for sub-table 1 (1D) + // - chi_part_tdndt_max : maximum chi for sub-table 1 (1D) + // - chi_part_tdndt_how_many : how many points to use for sub-table 1 (1D) + // - chi_part_tem_min : minimun chi for sub-table 2 (1D) + // - chi_part_tem_max : maximum chi for sub-table 2 (1D) + // - chi_part_tem_how_many : how many points to use for chi for sub-table 2 (2D) + // - prob_tem_how_many : how many points to use for the second axis of sub-table 2 (2D) + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real> ctrl; + + //Lookup table data + //---sub-table 1 (1D) + amrex::Gpu::ManagedDeviceVector<amrex::Real> KKfunc_coords; + amrex::Gpu::ManagedDeviceVector<amrex::Real> KKfunc_data; + //--- + + //---sub-table 2 (2D) + 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; + //______ +}; +//========================================================== + +#endif //WARPX_quantum_sync_engine_innards_h_ diff --git a/Source/QED/QuantumSyncEngineTableBuilder.H b/Source/QED/QuantumSyncEngineTableBuilder.H new file mode 100644 index 000000000..e70f5d02f --- /dev/null +++ b/Source/QED/QuantumSyncEngineTableBuilder.H @@ -0,0 +1,26 @@ +#ifndef WARPX_quantum_sync_engine_table_builder_h_ +#define WARPX_quantum_sync_engine_table_builder_h_ + +#include "QedWrapperCommons.H" +#include "QuantumSyncEngineInnards.H" + +//This includes only the definition of a simple datastructure +//used to control the Quantum Synchrotron engine. +#include <quantum_sync_engine_ctrl.h> + +/** + * A class which computes the lookup tables for the Quantum Synchrotron engine. + */ +class QuantumSynchrotronEngineTableBuilder{ +public: + /** + * Computes the tables. + * @param[in] ctrl control parameters to generate the tables + * @param[out] innards structure holding both a copy of ctrl and lookup tables data + */ + void compute_table + (picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real> ctrl, + QuantumSynchrotronEngineInnards& innards) const; +}; + +#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..51c3720f2 --- /dev/null +++ b/Source/QED/QuantumSyncEngineTableBuilder.cpp @@ -0,0 +1,58 @@ +#include "QuantumSyncEngineTableBuilder.H" + +//Include the full Quantum Synchrotron engine with table generation support +//(after some consistency tests). This requires to have a recent version +// of the Boost library. +#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> +//_______________________________________________ + +//Some handy aliases +using PicsarQuantumSynchrotronEngine = picsar::multi_physics:: + quantum_synchrotron_engine<amrex::Real, QedUtils::DummyStruct>; + +using PicsarQuantumSynchrotronCtrl = + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real>; +//_______________________________________________ + +void +QuantumSynchrotronEngineTableBuilder::compute_table + (PicsarQuantumSynchrotronCtrl ctrl, + QuantumSynchrotronEngineInnards& innards) const +{ + PicsarQuantumSynchrotronEngine qs_engine( + std::move(QedUtils::DummyStruct()), 1.0, ctrl); + + qs_engine.compute_dN_dt_lookup_table(); + qs_engine.compute_cumulative_phot_em_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); + innards.cum_distrib_coords_1.assign( + qs_innards_picsar.cum_distrib_table_coords_1_ptr, + qs_innards_picsar.cum_distrib_table_coords_1_ptr + + qs_innards_picsar.cum_distrib_table_coords_1_how_many); + innards.cum_distrib_coords_2.assign( + qs_innards_picsar.cum_distrib_table_coords_2_ptr, + qs_innards_picsar.cum_distrib_table_coords_2_ptr + + qs_innards_picsar.cum_distrib_table_coords_2_how_many); + innards.cum_distrib_data.assign( + qs_innards_picsar.cum_distrib_table_data_ptr, + qs_innards_picsar.cum_distrib_table_data_ptr + + qs_innards_picsar.cum_distrib_table_data_how_many); + //____ +} diff --git a/Source/QED/QuantumSyncEngineWrapper.H b/Source/QED/QuantumSyncEngineWrapper.H index a285dc8b6..1a6ffe4f3 100644 --- a/Source/QED/QuantumSyncEngineWrapper.H +++ b/Source/QED/QuantumSyncEngineWrapper.H @@ -2,18 +2,45 @@ #define WARPX_quantum_sync_engine_wrapper_h_ #include "QedWrapperCommons.H" +#include "QuantumSyncEngineInnards.H" -//QS ENGINE from PICSAR +#include <AMReX_Array.H> +#include <AMReX_Vector.H> +#include <AMReX_Gpu.H> + +//#define PXRMP_CORE_ONLY allows importing only the 'core functions' of the +//Quantum Synchrotron engine of the QED PICSAR library. #define PXRMP_CORE_ONLY -#include "quantum_sync_engine.hpp" +#include <quantum_sync_engine.hpp> + +//Lookup table building function is in a dedicated (optional) class to +//avoid including heavy dependencies if they are not needed. +#ifdef WARPX_QED_TABLE_GEN + #include "QuantumSyncEngineTableBuilder.H" +#endif + +#include <string> + +//Some handy aliases + +// The engine has two templated arguments: the numerical type +// and a random number generator. However, random numbers are not +// used to generate the lookup tables and the static member +// functions which are called in the functors do not use +// random numbers as well. Therefore, an empty "DummyStruct" +// can be passed as a template parameter. +using PicsarQuantumSynchrotronEngine = picsar::multi_physics:: + quantum_synchrotron_engine<amrex::Real, QedUtils::DummyStruct>; -using WarpXQuantumSynchrotronWrapper = - picsar::multi_physics::quantum_synchrotron_engine<amrex::Real, DummyStruct>; +using PicsarQuantumSynchrotronCtrl = + picsar::multi_physics::quantum_synchrotron_engine_ctrl<amrex::Real>; +//__________ // Functors ================================== -// These functors provide the core elementary functions of the library -// Can be included in GPU kernels +// 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. /** * Functor to initialize the optical depth of leptons for the @@ -22,33 +49,252 @@ using WarpXQuantumSynchrotronWrapper = class QuantumSynchrotronGetOpticalDepth { public: + /** + * Constructor does nothing because optical depth initialization + * does not require control parameters or lookup tables. + */ QuantumSynchrotronGetOpticalDepth () {}; + /** + * () operator is just a thin wrapper around a very simple function to + * generate the optical depth. It can be used on GPU. + */ AMREX_GPU_DEVICE AMREX_FORCE_INLINE - amrex::Real operator() () const + amrex::Real operator() () const noexcept { - return WarpXQuantumSynchrotronWrapper:: + //A random number in [0,1) should be provided as an argument. + return PicsarQuantumSynchrotronEngine:: internal_get_optical_depth(amrex::Random()); } }; //____________________________________________ +/** + * Functor to evolve the optical depth of leptons due to the + * Quantum Synchrotron process + */ +class QuantumSynchrotronEvolveOpticalDepth +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + QuantumSynchrotronEvolveOpticalDepth( + QuantumSynchrotronEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_KKfunc_size{r_innards.KKfunc_coords.size()}, + m_p_KKfunc_coords{r_innards.KKfunc_coords.dataPtr()}, + m_p_KKfunc_data{r_innards.KKfunc_data.dataPtr()} + {}; + + /** + * Evolves the optical depth. It can be used on GPU. + * @param[in] px,py,pz momentum components of the lepton (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] dt timestep (SI units) + * @param[in,out] opt_depth optical depth of the lepton. It is modified by the method. + * @return a flag which is 1 if optical depth becomes negative (i.e. a photon has to be generated). + */ + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + int 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 noexcept + { + bool has_event_happened{false}; + + //the library provides the time (< dt) at which the event occurs, but this + //feature won't be used in WarpX for now. + amrex::Real unused_event_time{0.0}; + + PicsarQuantumSynchrotronEngine:: + internal_evolve_opt_depth_and_determine_event( + px, py, pz, + ex, ey, ez, + bx, by, bz, + dt, opt_depth, + has_event_happened, unused_event_time, + m_dummy_lambda, + picsar::multi_physics::lookup_1d<amrex::Real>{ + m_KKfunc_size, + m_p_KKfunc_coords, + m_p_KKfunc_data}, + m_ctrl); + + return has_event_happened; + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarQuantumSynchrotronCtrl m_ctrl; + + //lookup table data + size_t m_KKfunc_size; + amrex::Real* m_p_KKfunc_coords; + amrex::Real* m_p_KKfunc_data; +}; + +/** + * Functor to generate a photon via the Quantum Synchrotron process + * and to update momentum accordingly + */ +class QuantumSynchrotronGeneratePhotonAndUpdateMomentum +{ +public: + /** + * Constructor acquires a reference to control parameters and + * lookup tables data. + * lookup_table uses non-owning vectors under the hood. So no new data + * allocations should be triggered on GPU + */ + QuantumSynchrotronGeneratePhotonAndUpdateMomentum( + QuantumSynchrotronEngineInnards& r_innards): + m_ctrl{r_innards.ctrl}, + m_cum_distrib_coords_1_size{r_innards.cum_distrib_coords_1.size()}, + m_cum_distrib_coords_2_size{r_innards.cum_distrib_coords_2.size()}, + m_p_distrib_coords_1{r_innards.cum_distrib_coords_1.data()}, + m_p_distrib_coords_2{r_innards.cum_distrib_coords_2.data()}, + m_p_cum_distrib_data{r_innards.cum_distrib_data.data()} + {}; + + /** + * Generates sampling (template parameter) photons according to Quantum Synchrotron process. + * It can be used on GPU. + * @param[in,out] px,py,pz momentum components of the lepton. They are modified (SI units) + * @param[in] ex,ey,ez electric field components (SI units) + * @param[in] bx,by,bz magnetic field components (SI units) + * @param[in] weight of the lepton (code units) + * @param[out] g_px,g_py,g_pz momenta of generated photons. Each array should have size=sampling (SI units) + * @param[out] g_weight weight of the generated photons. Array should have size=sampling (code units) + */ + template <size_t sampling> + AMREX_GPU_DEVICE + AMREX_FORCE_INLINE + void 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 noexcept + { + //[sampling] random numbers are needed + amrex::GpuArray<amrex::Real, sampling> + rand_zero_one_minus_epsi; + for(auto& el : rand_zero_one_minus_epsi) el = amrex::Random(); + + 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, + m_dummy_lambda, + picsar::multi_physics::lookup_2d<amrex::Real>{ + m_cum_distrib_coords_1_size, + m_p_distrib_coords_1, + m_cum_distrib_coords_2_size, + m_p_distrib_coords_2, + m_p_cum_distrib_data}, + m_ctrl, + rand_zero_one_minus_epsi.data()); + } + +private: + //laser wavelenght is not used with SI units + const amrex::Real m_dummy_lambda{1.0}; + + const PicsarQuantumSynchrotronCtrl m_ctrl; + + //lookup table data + size_t m_cum_distrib_coords_1_size; + size_t m_cum_distrib_coords_2_size; + amrex::Real* m_p_distrib_coords_1; + amrex::Real* m_p_distrib_coords_2; + amrex::Real* m_p_cum_distrib_data; +}; + // Factory class ============================= /** - * \brief Wrapper for the Quantum Synchrotron engine of the PICSAR library + * Wrapper for the Quantum Synchrotron engine of the PICSAR library */ class QuantumSynchrotronEngine { public: + /** + * Constructor requires no arguments. + */ QuantumSynchrotronEngine (); /** - * \brief Builds the functor to initialize the optical depth + * Builds the functor to initialize the optical depth */ QuantumSynchrotronGetOpticalDepth build_optical_depth_functor (); + + /** + * Builds the functor to evolve the optical depth + */ + QuantumSynchrotronEvolveOpticalDepth build_evolve_functor (); + + /** + * Builds the functor to generate photons + */ + QuantumSynchrotronGeneratePhotonAndUpdateMomentum build_phot_em_functor (); + + /** + * Checks if the optical tables are properly initialized + */ + bool are_lookup_tables_initialized () const; + + /** + * Init lookup tables from raw binary data. + * @param[in] raw_data a Vector of char + * @return true if it succeeds, false if it cannot parse raw_data + */ + bool init_lookup_tables_from_raw_data (const amrex::Vector<char>& raw_data); + + /** + * Export lookup tables data into a raw binary Vector + * @return the data in binary format. The Vector is empty if tables were + * not previously initialized. + */ + amrex::Vector<char> export_lookup_tables_data () const; + + /** + * Computes the lookup tables. It does nothing unless WarpX is compiled with QED_TABLE_GEN=TRUE + * @param[in] ctrl control params to generate the tables + */ + void compute_lookup_tables (PicsarQuantumSynchrotronCtrl ctrl); + + /** + * gets default (reasonable) values for the control parameters + * @return default control params to generate the tables + */ + PicsarQuantumSynchrotronCtrl get_default_ctrl() const; + +private: + bool m_lookup_tables_initialized = false; + + QuantumSynchrotronEngineInnards m_innards; + +//Table builing is available only if the libray is compiled with QED_TABLE_GEN=TRUE +#ifdef WARPX_QED_TABLE_GEN + QuantumSynchrotronEngineTableBuilder m_table_builder; +#endif + }; //============================================ diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp index b55149187..b2630dc4d 100644 --- a/Source/QED/QuantumSyncEngineWrapper.cpp +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -1,9 +1,16 @@ #include "QuantumSyncEngineWrapper.H" + +#include "QedTableParserHelperFunctions.H" + +#include <utility> + +using namespace std; +using namespace QedUtils; +using namespace amrex; + //This file provides a wrapper aroud the quantum_sync engine //provided by the PICSAR library -using namespace picsar::multi_physics; - // Factory class ============================= QuantumSynchrotronEngine::QuantumSynchrotronEngine (){} @@ -14,4 +21,167 @@ QuantumSynchrotronEngine::build_optical_depth_functor () return QuantumSynchrotronGetOpticalDepth(); } +QuantumSynchrotronEvolveOpticalDepth QuantumSynchrotronEngine::build_evolve_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return QuantumSynchrotronEvolveOpticalDepth(m_innards); +} + +QuantumSynchrotronGeneratePhotonAndUpdateMomentum QuantumSynchrotronEngine::build_phot_em_functor () +{ + AMREX_ALWAYS_ASSERT(m_lookup_tables_initialized); + + return QuantumSynchrotronGeneratePhotonAndUpdateMomentum(m_innards); + +} + +bool QuantumSynchrotronEngine::are_lookup_tables_initialized () const +{ + return m_lookup_tables_initialized; +} + +bool +QuantumSynchrotronEngine::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, m_innards.ctrl.chi_part_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tdndt_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tdndt_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tdndt_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tdndt_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tdndt_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tdndt_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tem_min, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tem_min)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tem_max, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tem_max)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.chi_part_tem_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.chi_part_tem_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + tie(is_ok, m_innards.ctrl.prob_tem_how_many, p_data) = + parse_raw_data<decltype(m_innards.ctrl.prob_tem_how_many)>( + p_data, p_last); + if(!is_ok) return false; + + //___________________________ + + //Data + Vector<Real> tndt_coords(m_innards.ctrl.chi_part_tdndt_how_many); + Vector<Real> tndt_data(m_innards.ctrl.chi_part_tdndt_how_many); + Vector<Real> cum_tab_coords1(m_innards.ctrl.chi_part_tem_how_many); + Vector<Real> cum_tab_coords2(m_innards.ctrl.prob_tem_how_many); + Vector<Real> cum_tab_data(m_innards.ctrl.chi_part_tem_how_many* + m_innards.ctrl.prob_tem_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; + m_innards.KKfunc_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; + m_innards.KKfunc_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; + m_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; + m_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; + m_innards.cum_distrib_data.assign( + cum_tab_data.begin(), cum_tab_data.end()); + + //___________________________ + m_lookup_tables_initialized = true; + + return true; +} + +Vector<char> QuantumSynchrotronEngine::export_lookup_tables_data () const +{ + Vector<char> res{}; + + if(!m_lookup_tables_initialized) + return res; + + add_data_to_vector_char(&m_innards.ctrl.chi_part_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tdndt_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tdndt_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tdndt_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tem_min, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tem_max, 1, res); + add_data_to_vector_char(&m_innards.ctrl.chi_part_tem_how_many, 1, res); + add_data_to_vector_char(&m_innards.ctrl.prob_tem_how_many, 1, res); + + add_data_to_vector_char(m_innards.KKfunc_coords.data(), + m_innards.KKfunc_coords.size(), res); + add_data_to_vector_char(m_innards.KKfunc_data.data(), + m_innards.KKfunc_data.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_1.data(), + m_innards.cum_distrib_coords_1.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_coords_2.data(), + m_innards.cum_distrib_coords_2.size(), res); + add_data_to_vector_char(m_innards.cum_distrib_data.data(), + m_innards.cum_distrib_data.size(), res); + + return res; +} + +PicsarQuantumSynchrotronCtrl +QuantumSynchrotronEngine::get_default_ctrl() const +{ + return PicsarQuantumSynchrotronCtrl(); +} + +void QuantumSynchrotronEngine::compute_lookup_tables ( + PicsarQuantumSynchrotronCtrl ctrl) +{ +#ifdef WARPX_QED_TABLE_GEN + m_table_builder.compute_table(ctrl, m_innards); + m_lookup_tables_initialized = true; +#endif +} + //============================================ diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index fd6e72dc6..ab28c5446 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -2,6 +2,8 @@ #include <AMReX_Vector.H> #include <AMReX_MultiFab.H> +#include <string> + void ReadBoostedFrameParameters(amrex::Real& gamma_boost, amrex::Real& beta_boost, amrex::Vector<int>& boost_direction); @@ -9,3 +11,13 @@ void ConvertLabParamsToBoost(); void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax); + +namespace WarpXUtilIO{ + /** + * A helper function to write binary data on disk. + * @param[in] filename where to write + * @param[in] data Vector containing binary data to write on disk + * return true if it succeeds, false otherwise + */ + bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data); +} diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 4b11eb69d..8764a09c6 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -1,10 +1,11 @@ -#include <cmath> - #include <WarpXUtil.H> #include <WarpXConst.H> #include <AMReX_ParmParse.H> #include <WarpX.H> +#include <cmath> +#include <fstream> + using namespace amrex; void ReadBoostedFrameParameters(Real& gamma_boost, Real& beta_boost, @@ -152,3 +153,14 @@ void NullifyMF(amrex::MultiFab& mf, int lev, amrex::Real zmin, amrex::Real zmax) } } } + +namespace WarpXUtilIO{ + bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data) + { + std::ofstream of{filename, std::ios::binary}; + of.write(data.data(), data.size()); + of.close(); + return of.good(); + } +} + |