diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 119dfab5d99..e57abea36ce 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -33,6 +33,7 @@ #include "Utils/WarpXUtil.H" #include "Utils/WarpXConst.H" #include "Utils/WarpXProfilerWrapper.H" +#include "Particles/Radiation/RadiationHandler.H" #include #include @@ -245,6 +246,8 @@ WarpX::Evolve (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1], *Bfield_aux[lev][2]); } + // For now only on lower level + mypc->Dump_radiations(dt[0]); is_synchronized = true; } @@ -426,8 +429,20 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("particlescraper"); ExecutePythonCallback("beforedeposition"); + //Save particles' old momenta + mypc->RecordOldMomenta(); + PushParticlesandDeposit(cur_time); + //Radiation contribution at each timestep + //Only level 0 is supported + mypc->doRadiation(dt[0],cur_time); + + + + + + ExecutePythonCallback("afterdeposition"); // Synchronize J and rho: diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp index 640cf1b82af..851d92ff29a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmComoving.cpp @@ -3,7 +3,8 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" + +#include #include #include @@ -22,6 +23,7 @@ #if WARPX_USE_PSATD using namespace amrex; +using namespace ablastr::math; PsatdAlgorithmComoving::PsatdAlgorithmComoving (const SpectralKSpace& spectral_kspace, const DistributionMapping& dm, @@ -75,7 +77,7 @@ PsatdAlgorithmComoving::pushSpectralFields (SpectralFieldData& f) const const amrex::Box& bx = f.fields[mfi].box(); // Extract arrays for the fields to be updated - const amrex::Array4 fields = f.fields[mfi].array(); + const amrex::Array4 fields = f.fields[mfi].array(); // Extract arrays for the coefficients const amrex::Array4 C_arr = C_coef [mfi].array(); diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp index 6d4d8613940..36f5c3c99b6 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp @@ -9,7 +9,8 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" + +#include #include #include @@ -27,6 +28,7 @@ #if WARPX_USE_PSATD using namespace amrex::literals; +using namespace ablastr::math; PsatdAlgorithmFirstOrder::PsatdAlgorithmFirstOrder( const SpectralKSpace& spectral_kspace, diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp index 167d7565bc0..c00ab2ff79e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp @@ -9,7 +9,7 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" +#include #include #include @@ -28,6 +28,7 @@ #if WARPX_USE_PSATD using namespace amrex; +using namespace ablastr::math; PsatdAlgorithmJConstantInTime::PsatdAlgorithmJConstantInTime( const SpectralKSpace& spectral_kspace, diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp index 92e35c8c03b..f9280798fd0 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJLinearInTime.cpp @@ -8,7 +8,8 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" + +#include #include #include @@ -27,6 +28,7 @@ #if WARPX_USE_PSATD using namespace amrex::literals; +using namespace ablastr::math; PsatdAlgorithmJLinearInTime::PsatdAlgorithmJLinearInTime( const SpectralKSpace& spectral_kspace, diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp index 92b8b3e900d..715b29fd62c 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmPml.cpp @@ -11,7 +11,8 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" + +#include #include #include @@ -29,6 +30,7 @@ #if WARPX_USE_PSATD using namespace amrex; +using namespace ablastr::math; PsatdAlgorithmPml::PsatdAlgorithmPml( const SpectralKSpace& spectral_kspace, diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H index 4af24cfa559..3cae6502469 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H @@ -8,11 +8,12 @@ #define WARPX_SPECTRAL_BASE_ALGORITHM_H_ #include "FieldSolver/SpectralSolver/SpectralKSpace.H" -#include "Utils/WarpX_Complex.H" #include "FieldSolver/SpectralSolver/SpectralFieldData_fwd.H" #include "FieldSolver/SpectralSolver/SpectralFieldData.H" +#include + #include #include #include @@ -80,7 +81,7 @@ class SpectralBaseAlgorithm using SpectralRealCoefficients = \ amrex::FabArray< amrex::BaseFab >; using SpectralComplexCoefficients = \ - amrex::FabArray< amrex::BaseFab >; + amrex::FabArray< amrex::BaseFab >; /** * \brief Constructor diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp index d635a5debe3..99a98900c9e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp @@ -7,7 +7,8 @@ #include "SpectralBaseAlgorithm.H" #include "FieldSolver/SpectralSolver/SpectralFieldData.H" -#include "Utils/WarpX_Complex.H" + +#include #include #include @@ -23,6 +24,7 @@ #include using namespace amrex; +using namespace ablastr::math; /** * \brief Constructor diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index 8ced43ced94..065fbaecc4a 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -11,8 +11,8 @@ #include "SpectralFieldData_fwd.H" #include "SpectralKSpace.H" -#include "Utils/WarpX_Complex.H" +#include #include #include @@ -28,7 +28,7 @@ #include // Declare type for spectral fields -using SpectralField = amrex::FabArray< amrex::BaseFab >; +using SpectralField = amrex::FabArray< amrex::BaseFab >; class SpectralFieldIndex { diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 35d37df6038..072f32dab5d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -11,6 +11,8 @@ #include "Utils/WarpXUtil.H" #include "WarpX.H" +#include + #include #include #include @@ -32,6 +34,7 @@ #if WARPX_USE_PSATD using namespace amrex; +using namespace ablastr::math; SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H index 2a4dd3d1580..101f7909fe1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.H +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.H @@ -10,7 +10,7 @@ #include "SpectralKSpace_fwd.H" -#include "Utils/WarpX_Complex.H" +#include #include #include @@ -29,7 +29,7 @@ using RealKVector = amrex::Gpu::DeviceVector; using KVectorComponent = amrex::LayoutData< RealKVector >; using SpectralShiftFactor = amrex::LayoutData< - amrex::Gpu::DeviceVector >; + amrex::Gpu::DeviceVector >; // Indicate the type of correction "shift" factor to apply // when the FFT is performed from/to a cell-centered grid in real space. diff --git a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp index 91fe2953668..42078b2d8cd 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp @@ -11,6 +11,8 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" +#include + #include #include #include @@ -27,6 +29,7 @@ #include using namespace amrex; +using namespace ablastr::math; /* \brief Initialize k space object. * diff --git a/Source/Laser/LaserProfiles.H b/Source/Laser/LaserProfiles.H index 3aea88ac956..61eba2edb5f 100644 --- a/Source/Laser/LaserProfiles.H +++ b/Source/Laser/LaserProfiles.H @@ -7,6 +7,8 @@ #ifndef WARPX_LaserProfiles_H_ #define WARPX_LaserProfiles_H_ +#include + #include #include #include @@ -22,8 +24,6 @@ #include #include -#include "Utils/WarpX_Complex.H" - namespace WarpXLaserProfiles { /** Common laser profile parameters @@ -371,7 +371,7 @@ private: /** Index of the last timestep in memory */ int last_time_index; /** lasy field data */ - amrex::Gpu::DeviceVector E_lasy_data; + amrex::Gpu::DeviceVector E_lasy_data; /** binary field data */ amrex::Gpu::DeviceVector E_binary_data; /** This parameter is subtracted to simulation time before interpolating field data in file (either lasy or binary). diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp index 876a56d537e..95e051b6a2d 100644 --- a/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp +++ b/Source/Laser/LaserProfilesImpl/LaserProfileFieldFunction.cpp @@ -8,7 +8,6 @@ #include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" -#include "Utils/WarpX_Complex.H" #include #include diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp index 0fdb45c64f8..3787ec43118 100644 --- a/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp +++ b/Source/Laser/LaserProfilesImpl/LaserProfileFromFile.cpp @@ -9,9 +9,9 @@ #include "Utils/Algorithms/LinearInterpolation.H" #include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" -#include "Utils/WarpX_Complex.H" #include "Utils/WarpXConst.H" +#include #include #include @@ -48,6 +48,7 @@ #endif using namespace amrex; +using namespace ablastr::math; void WarpXLaserProfiles::FromFileLaserProfile::init ( diff --git a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp index 96d61d920f1..e2db520c014 100644 --- a/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp +++ b/Source/Laser/LaserProfilesImpl/LaserProfileGaussian.cpp @@ -10,12 +10,12 @@ #include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" -#include "Utils/WarpX_Complex.H" + +#include #include #include #include -#include #include #include #include @@ -28,6 +28,7 @@ #include using namespace amrex; +using namespace ablastr::math; void WarpXLaserProfiles::GaussianLaserProfile::init ( diff --git a/Source/Particles/CMakeLists.txt b/Source/Particles/CMakeLists.txt index 67af14ef889..01800fa7cc8 100644 --- a/Source/Particles/CMakeLists.txt +++ b/Source/Particles/CMakeLists.txt @@ -21,5 +21,6 @@ add_subdirectory(ElementaryProcess) add_subdirectory(Gather) add_subdirectory(ParticleCreation) #add_subdirectory(Pusher) +add_subdirectory(Radiation) add_subdirectory(Resampling) add_subdirectory(Sorting) diff --git a/Source/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index d0db678dfda..128a0edabc4 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -13,8 +13,9 @@ #include "Particles/Pusher/GetAndSetPosition.H" #include "Particles/ShapeFactors.H" #include "Utils/WarpXAlgorithmSelection.H" + #ifdef WARPX_DIM_RZ -# include "Utils/WarpX_Complex.H" +# include #endif #include @@ -52,6 +53,10 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, { using namespace amrex; +#ifdef WARPX_DIM_RZ + using namespace ablastr::math; +#endif + #if !defined(AMREX_USE_GPU) amrex::ignore_unused(cost, load_balance_costs_update_algo); #endif @@ -260,6 +265,9 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio ) { using namespace amrex; +#ifdef WARPX_DIM_RZ + using namespace ablastr::math; +#endif const auto *permutation = a_bins.permutationPtr(); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index ac377a6799d..c42708fcf0d 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -15,12 +15,13 @@ #include "Utils/TextMsg.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" -#ifdef WARPX_DIM_RZ -# include "Utils/WarpX_Complex.H" -#endif #include "WarpX.H" // todo: remove include and pass globals as args +#ifdef WARPX_DIM_RZ +# include +#endif + #include #include #include @@ -71,6 +72,9 @@ void doDepositionShapeN (const GetParticlePosition& GetPosition, long load_balance_costs_update_algo) { using namespace amrex::literals; +#ifdef WARPX_DIM_RZ + using namespace ablastr::math; +#endif #if !defined(WARPX_DIM_RZ) amrex::ignore_unused(n_rz_azimuthal_modes); @@ -599,6 +603,9 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, { using namespace amrex; using namespace amrex::literals; +#ifdef WARPX_DIM_RZ + using namespace ablastr::math; +#endif #if !defined(WARPX_DIM_RZ) ignore_unused(n_rz_azimuthal_modes); diff --git a/Source/Particles/Deposition/SharedDepositionUtils.H b/Source/Particles/Deposition/SharedDepositionUtils.H index e28835a57df..3b46c340f35 100644 --- a/Source/Particles/Deposition/SharedDepositionUtils.H +++ b/Source/Particles/Deposition/SharedDepositionUtils.H @@ -11,8 +11,9 @@ #include "Particles/ShapeFactors.H" #include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpXConst.H" + #ifdef WARPX_DIM_RZ -# include "Utils/WarpX_Complex.H" +# include #endif #include @@ -86,6 +87,8 @@ void depositComponent (const GetParticlePosition& GetPosition, #if !defined(WARPX_DIM_RZ) amrex::ignore_unused(n_rz_azimuthal_modes); +#else + using namespace ablastr::math; #endif // Whether ion_lev is a null pointer (do_ionization=0) or a real pointer diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 1b5cede3c6f..be17670ffea 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -11,7 +11,8 @@ #include "Particles/Gather/GetExternalFields.H" #include "Particles/Pusher/GetAndSetPosition.H" #include "Particles/ShapeFactors.H" -#include "Utils/WarpX_Complex.H" + +#include #include @@ -62,6 +63,7 @@ void doGatherShapeN (const amrex::ParticleReal xp, const int n_rz_azimuthal_modes) { using namespace amrex; + using namespace ablastr::math; #if defined(WARPX_DIM_XZ) amrex::ignore_unused(yp); diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 58cbe11a980..e873088b7c9 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -16,6 +16,7 @@ include $(WARPX_HOME)/Source/Particles/Sorting/Make.package include $(WARPX_HOME)/Source/Particles/ParticleCreation/Make.package include $(WARPX_HOME)/Source/Particles/ElementaryProcess/Make.package include $(WARPX_HOME)/Source/Particles/Collision/Make.package +include $(WARPX_HOME)/Source/Particles/Radiation/Make.package include $(WARPX_HOME)/Source/Particles/Filter/Make.package include $(WARPX_HOME)/Source/Particles/Resampling/Make.package diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index c32ff07464d..791f7c77329 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -20,6 +20,7 @@ # include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper_fwd.H" #endif #include "PhysicalParticleContainer.H" +#include "Particles/Radiation/RadiationHandler_fwd.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" #include "WarpXParticleContainer.H" @@ -68,7 +69,7 @@ public: MultiParticleContainer (amrex::AmrCore* amr_core); - ~MultiParticleContainer() = default; + ~MultiParticleContainer(); MultiParticleContainer (MultiParticleContainer const &) = delete; MultiParticleContainer& operator= (MultiParticleContainer const & ) = delete; @@ -329,6 +330,11 @@ public: const amrex::MultiFab& Bz); #endif + void doRadiation(amrex::Real dt, amrex::Real cur_time); + + void RecordOldMomenta(); + + void Dump_radiations(amrex::Real dt); [[nodiscard]] int getSpeciesID (std::string product_str) const; protected: @@ -480,6 +486,12 @@ protected: amrex::Real m_qed_schwinger_zmin = std::numeric_limits::lowest(); amrex::Real m_qed_schwinger_zmax = std::numeric_limits::max(); + /** + * For Radiation module + */ + bool m_at_least_one_particle_radiate = 0; + int m_nspecies_radiate; + #endif private: @@ -527,6 +539,8 @@ private: void CheckQEDProductSpecies(); #endif - +// Radiation Handler +std::unique_ptr m_p_radiation_handler; +bool m_at_least_one_has_radiation = false; }; #endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 213a5344ee8..7970511dbfd 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -29,6 +29,7 @@ #include "Particles/ParticleCreation/SmartUtils.H" #include "Particles/PhotonParticleContainer.H" #include "Particles/PhysicalParticleContainer.H" +#include "Particles/Radiation/RadiationHandler.H" #include "Particles/RigidInjectedParticleContainer.H" #include "Particles/WarpXParticleContainer.H" #include "SpeciesPhysicalProperties.H" @@ -122,9 +123,24 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) // Setup particle collisions collisionhandler = std::make_unique(this); - + //Initialization of the radiation + for (auto& s : allcontainers) { + if (s->has_radiation()) { + m_at_least_one_has_radiation = true; + + const auto& lev_0_geom = amr_core->Geom(0); + auto center = amrex::Array{}; + for(int idim=0; idim(center); + break; + } + } } +MultiParticleContainer::~MultiParticleContainer() {} + void MultiParticleContainer::ReadParameters () { @@ -367,8 +383,7 @@ MultiParticleContainer::ReadParameters () m_qed_schwinger_threshold_poisson_gaussian); utils::parser::queryWithParser( pp_qed_schwinger, "xmin", m_qed_schwinger_xmin); - utils::parser::queryWithParser( - pp_qed_schwinger, "xmax", m_qed_schwinger_xmax); + #if defined(WARPX_DIM_3D) utils::parser::queryWithParser( pp_qed_schwinger, "ymin", m_qed_schwinger_ymin); @@ -384,7 +399,6 @@ MultiParticleContainer::ReadParameters () initialized = true; } } - WarpXParticleContainer& MultiParticleContainer::GetParticleContainerFromName (const std::string& name) const { @@ -915,6 +929,14 @@ MultiParticleContainer::doFieldIonization (int lev, } } +void MultiParticleContainer::RecordOldMomenta(){ + for (auto& pc : allcontainers) { + if (pc->has_radiation()){ + pc->RecordOldMomenta(); + } + } +} + void MultiParticleContainer::doCollisions ( Real cur_time, amrex::Real dt ) { @@ -955,6 +977,23 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vectorhas_radiation()){ + m_p_radiation_handler->add_radiation_contribution(dt,pc,cur_time); + } + } + } +} + +void MultiParticleContainer::Dump_radiations(const amrex::Real dt){ + if (m_at_least_one_has_radiation){ + m_p_radiation_handler->Integral_overtime(dt); + m_p_radiation_handler->gather_and_write_radiation("Radiation"); + } +} #ifdef WARPX_QED void MultiParticleContainer::InitQED () { diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index 3be0886425d..394b808ad8b 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -177,6 +177,23 @@ public: /** Return the name-to-index map for the runtime-time integer components */ [[nodiscard]] std::map getParticleRuntimeiComps () const noexcept { return particle_runtime_icomps;} + /** Get the index of a real component + * + * @param name Name of the new component + */ + int GetRealCompIndex (const std::string& name) + { + return particle_comps.at(name); + } + /** Initialize value at a component + * + * @param name Name of the new component + */ + // void InitializeComp (const std::string& name) + //{ + // particle_runtime_comps.at(name); + //} + protected: std::map particle_comps; std::map particle_icomps; diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 04c9682486c..0453a626071 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -358,6 +358,8 @@ public: PairGenerationFilterFunc getPairGenerationFilterFunc (); #endif + bool has_radiation() const override {return m_do_radiation;} + protected: std::string species_name; std::vector> plasma_injectors; @@ -386,6 +388,7 @@ protected: //When true PhysicalParticleContainer tries to use a pusher including //radiation reaction bool do_classical_radiation_reaction = false; + bool m_do_radiation = false; // A flag to enable saving of the previous timestep positions bool m_save_previous_position = false; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 729abd6f092..645f90f0937 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -347,6 +347,11 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //check if Radiation Reaction is enabled and do consistency checks pp_species_name.query("do_classical_radiation_reaction", do_classical_radiation_reaction); + + //check if Radiations are enables for species and add old momentum attribute + pp_species_name.query("do_radiation", m_do_radiation); + if (m_do_radiation==1) {AllocateOldMomenta();} + //if the species is not a lepton, do_classical_radiation_reaction //should be false WARPX_ALWAYS_ASSERT_WITH_MESSAGE( diff --git a/Source/Particles/Radiation/CMakeLists.txt b/Source/Particles/Radiation/CMakeLists.txt new file mode 100644 index 00000000000..b7b0d37d608 --- /dev/null +++ b/Source/Particles/Radiation/CMakeLists.txt @@ -0,0 +1,7 @@ +foreach(D IN LISTS WarpX_DIMS) + warpx_set_suffix_dims(SD ${D}) + target_sources(lib_${SD} + PRIVATE + RadiationHandler.cpp + ) +endforeach() diff --git a/Source/Particles/Radiation/Make.package b/Source/Particles/Radiation/Make.package new file mode 100644 index 00000000000..1c3a1020746 --- /dev/null +++ b/Source/Particles/Radiation/Make.package @@ -0,0 +1,3 @@ +CEXE_sources += RadiationHandler.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Radiation diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H new file mode 100644 index 00000000000..1b1a2c8f3a0 --- /dev/null +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -0,0 +1,66 @@ + +/* Copyright 2023 Thomas Clark, Luca Fedeli + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_PARTICLES_RADIATION_H +#define WARPX_PARTICLES_RADIATION_H + +#include "RadiationHandler_fwd.H" + +#include "Particles/WarpXParticleContainer_fwd.H" + +#include + +#include +#include +#include + +#include +#include + +/* \brief CollisionHandler is a class that contains the + * calculation of radiation for particles at each frequencies and angles. + */ +class RadiationHandler +{ +public: + RadiationHandler (const amrex::Array& center); + + /* Perform the calculation of the radiation */ + //void doRadiation (amrex::Real dt, MultiParticleContainer* mypc); + void add_radiation_contribution(const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time); + + void gather_and_write_radiation(const std::string& filename); + + void Integral_overtime(const amrex::Real dt); + +private: + + // Frequency range + amrex::Array m_omega_range; + int m_omega_points; + amrex::Gpu::DeviceVector m_omegas; + + // Dimensions of the detector + amrex::Array m_theta_range; + + amrex::Real m_det_distance; + amrex::Array m_det_pts; + amrex::Array m_det_direction; + amrex::Array m_det_orientation; + + amrex::Gpu::DeviceVector det_pos_x; + amrex::Gpu::DeviceVector det_pos_y; + amrex::Gpu::DeviceVector det_pos_z; + amrex::Gpu::DeviceVector det_pos_theta; + amrex::Gpu::DeviceVector det_pos_phi; + + amrex::Gpu::DeviceVector m_radiation_data; + amrex::Gpu::DeviceVector m_radiation_calculation; + +}; +#endif // WARPX_PARTICLES_RADIATION_H diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp new file mode 100644 index 00000000000..c57047631c0 --- /dev/null +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -0,0 +1,445 @@ +/* Copyright 2023 Thomas Clark, Luca Fedeli + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "RadiationHandler.H" + +#include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/WarpXParticleContainer.H" +#include "Utils/Parser/ParserUtils.H" +#include "Utils/WarpXConst.H" +#include "Utils/TextMsg.H" + +#include + +#ifdef AMREX_USE_OMP +# include +#endif + +#include +#include + +using namespace amrex; +using namespace ablastr::math; +using namespace utils::parser; + +enum class RadiationType +{ + cartesian, + spherical +}; +auto const m_radiation_type = std::map{ + {"cartesian", RadiationType::cartesian}, + {"spherical", RadiationType::spherical} +}; + +namespace +{ + auto compute_detector_positions( + const amrex::Array& center, + const amrex::Array& direction, + const amrex::Real distance, + const amrex::Array& orientation, + const amrex::Array& det_points, + const amrex::Array& theta_range, + const std::string radiation_type) + { + + const auto how_many = det_points[0]*det_points[1]; + + auto host_det_x = amrex::Vector(how_many); + auto host_det_y = amrex::Vector(how_many); + auto host_det_z = amrex::Vector(how_many); + auto host_det_theta = amrex::Vector(how_many); + auto host_det_phi = amrex::Vector(how_many); + + if(m_radiation_type.at(radiation_type) == RadiationType::cartesian){ + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + direction[0]*orientation[0] + + direction[1]*orientation[1] + + direction[2]*orientation[2] == 0, + "Radiation detector orientation cannot be aligned with detector direction"); + + auto u = amrex::Array{ + orientation[0] - direction[0]*orientation[0], + orientation[1] - direction[1]*orientation[1], + orientation[2] - direction[2]*orientation[2]}; + const auto one_over_u = 1.0_rt/std::sqrt(u[0]*u[0]+u[1]*u[1]+u[2]*u[2]); + u[0] *= one_over_u; + u[1] *= one_over_u; + u[2] *= one_over_u; + + auto v = amrex::Array{ + direction[1]*u[2]-direction[2]*u[1], + direction[2]*u[0]-direction[0]*u[2], + direction[0]*u[1]-direction[1]*u[0]}; + const auto one_over_v = 1.0_rt/std::sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); + v[0] *= one_over_v; + v[1] *= one_over_v; + v[2] *= one_over_v; + + auto us = amrex::Vector(det_points[0]); + const auto ulim = distance*std::tan(theta_range[0]*0.5_rt); + amrex::linspace(us.begin(), us.end(), -ulim,ulim); + + auto vs = amrex::Vector(det_points[1]); + const auto vlim = distance*std::tan(theta_range[1]*0.5_rt); + amrex::linspace(vs.begin(), vs.end(), -vlim, vlim); + + const auto one_over_direction = 1.0_rt/std::sqrt( + direction[0]*direction[0]+direction[1]*direction[1]+direction[2]*direction[2]); + const auto norm_direction = amrex::Array{ + direction[0]*one_over_direction, + direction[1]*one_over_direction, + direction[2]*one_over_direction}; + + for (int i = 0; i < det_points[0]; ++i) + { + for (int j = 0; j < det_points[1]; ++j) + { + auto x = distance * norm_direction[0] + us[i]*u[0] + vs[j]*v[0]; + auto y = distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; + auto z = distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; + + host_det_x[i*det_points[1] + j] = center[0] + x; + host_det_y[i*det_points[1] + j] = center[1] + y; + host_det_z[i*det_points[1] + j] = center[2] + z; + + host_det_theta[i*det_points[1] + j] = std::acos(z/distance); + host_det_phi[i*det_points[1] + j] = std::atan2(y, x) + ablastr::constant::math::pi; + + } + } + } + + if(m_radiation_type.at(radiation_type) == RadiationType::spherical){ + + const auto Ntheta = det_points[0]; + const auto Nphi = det_points[1]; + const auto dtheta = theta_range[0]/Ntheta; + const auto dphi = theta_range[1]/Nphi; + const auto thetaOrigin = std::acos(direction[2]) - Ntheta/2*dtheta; + const auto phiOrigin = std::atan2(direction[1], direction[0]) - Nphi/2*dphi; + + for (int i = 0; i < det_points[0]; ++i) + { + for (int j = 0; j < det_points[1]; ++j) + { + auto theta = thetaOrigin + i*dtheta; + auto phi = phiOrigin + j*dphi; + + //theta is in [0, pi] and phi is in [0, 2*pi] + if(theta < 0){ + theta *= -1; + phi += ablastr::constant::math::pi; + } + phi -= 2*ablastr::constant::math::pi*std::floor(phi/(2*ablastr::constant::math::pi)); + + + host_det_x[i*det_points[1] + j] = center[0] + distance*std::sin(theta)*std::cos(phi); + host_det_y[i*det_points[1] + j] = center[1] + distance*std::sin(theta)*std::sin(phi); + host_det_z[i*det_points[1] + j] = center[2] + distance*std::cos(theta); + + host_det_theta[i*det_points[1] + j] = theta; + host_det_phi[i*det_points[1] + j] = phi; + } + } + } + + auto gpu_det_x = amrex::Gpu::DeviceVector(how_many); + auto gpu_det_y = amrex::Gpu::DeviceVector(how_many); + auto gpu_det_z = amrex::Gpu::DeviceVector(how_many); + auto gpu_det_theta = amrex::Gpu::DeviceVector(how_many); + auto gpu_det_phi = amrex::Gpu::DeviceVector(how_many); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_det_x.begin(), host_det_x.end(), gpu_det_x.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_det_y.begin(), host_det_y.end(), gpu_det_y.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_det_z.begin(), host_det_z.end(), gpu_det_z.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_det_theta.begin(), host_det_theta.end(), gpu_det_theta.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + host_det_phi.begin(), host_det_phi.end(), gpu_det_phi.begin()); + + + amrex::Gpu::Device::streamSynchronize(); + + return std::make_tuple(gpu_det_x, gpu_det_y, gpu_det_z, gpu_det_theta, gpu_det_phi); + + } +} + + +RadiationHandler::RadiationHandler(const amrex::Array& center) +{ +#if defined WARPX_DIM_RZ + WARPX_ABORT_WITH_MESSAGE("Radiation is not supported yet with RZ."); +#endif + + // Read in radiation input + const amrex::ParmParse pp_radiation("radiation"); + + //type of detector + std::string radiation_type = "spherical"; + pp_radiation.query("detector_type", radiation_type); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_radiation_type.find(radiation_type) != m_radiation_type.end(), radiation_type + " detector type is not supported yet"); + + //Resolution in frequency of the detector + auto omega_range = std::vector(2); + getArrWithParser(pp_radiation, "omega_range", omega_range); + std::copy(omega_range.begin(), omega_range.end(), m_omega_range.begin()); + getWithParser(pp_radiation, "omega_points", m_omega_points); + + //Angle theta AND phi + auto theta_range = std::vector(2); + getArrWithParser(pp_radiation, "angle_aperture", theta_range); + std::copy(theta_range.begin(), theta_range.end(), m_theta_range.begin()); + + //Detector parameters + auto det_pts = std::vector(2); + getArrWithParser(pp_radiation, "detector_number_points", det_pts); + std::copy(det_pts.begin(), det_pts.end(), m_det_pts.begin()); + + auto det_direction = std::vector(3); + getArrWithParser(pp_radiation, "detector_direction", det_direction); + std::copy(det_direction.begin(), det_direction.end(), m_det_direction.begin()); + + auto det_orientation = std::vector(3); + getArrWithParser(pp_radiation, "detector_orientation", det_orientation); + std::copy(det_orientation.begin(), det_orientation.end(), m_det_orientation.begin()); + + getWithParser(pp_radiation, "detector_distance", m_det_distance); + + std::tie(det_pos_x, det_pos_y, det_pos_z, det_pos_theta, det_pos_phi) = compute_detector_positions( + center, m_det_direction, m_det_distance, + m_det_orientation, m_det_pts, m_theta_range, radiation_type); + + constexpr auto ncomp = 3; + m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); + + auto t_omegas = amrex::Vector(m_omega_points); + amrex::linspace(t_omegas.begin(), t_omegas.end(), m_omega_range[0], m_omega_range[1]); + m_omegas = amrex::Gpu::DeviceVector(m_omega_points); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + t_omegas.begin(), t_omegas.end(), m_omegas.begin()); + amrex::Gpu::Device::streamSynchronize(); +} + + +void RadiationHandler::add_radiation_contribution + (const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time) + { + for (int lev = 0; lev <= pc->finestLevel(); ++lev) + { +#ifdef AMREX_USE_OMP + #pragma omp parallel +#endif + { + + for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) + { + long const np = pti.numParticles(); + const auto& attribs = pti.GetAttribs(); + const auto* p_w = attribs[PIdx::w].data(); + const auto* p_ux = attribs[PIdx::ux].data(); + const auto* p_uy = attribs[PIdx::uy].data(); + const auto* p_uz = attribs[PIdx::uz].data(); + + const auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + auto& part = pc->GetParticles(lev)[index]; + auto& soa = part.GetStructOfArrays(); + + const auto* p_ux_old = soa.GetRealData(pc->GetRealCompIndex("old_u_x")).data(); + const auto* p_uy_old = soa.GetRealData(pc->GetRealCompIndex("old_u_y")).data(); + const auto* p_uz_old = soa.GetRealData(pc->GetRealCompIndex("old_u_z")).data(); + + auto GetPosition = GetParticlePosition(pti); + auto const q = pc->getCharge(); + + const auto how_many_det_pos = static_cast(det_pos_x.size()); + + const auto p_omegas = m_omegas.dataPtr(); + + const auto p_det_pos_x = det_pos_x.dataPtr(); + const auto p_det_pos_y = det_pos_y.dataPtr(); + const auto p_det_pos_z = det_pos_z.dataPtr(); + + auto p_radiation_data = m_radiation_data.dataPtr(); + + const auto& omega_points = m_omega_points; + + constexpr auto c = PhysConst::c; + constexpr auto inv_c = 1._prt/(PhysConst::c); + constexpr auto inv_c2 = 1._prt/(PhysConst::c*PhysConst::c); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip){ + amrex::ParticleReal xp, yp, zp; + GetPosition.AsStored(ip,xp, yp, zp); + + const auto ux = 0.5_prt*(p_ux[ip] + p_ux_old[ip]); + const auto uy = 0.5_prt*(p_uy[ip] + p_uy_old[ip]); + const auto uz = 0.5_prt*(p_uz[ip] + p_uz_old[ip]); + auto const u2 = ux*ux + uy*uy + uz*uz; + + auto const one_over_gamma = 1._prt/std::sqrt(1.0_rt + u2*inv_c2); + + auto const one_over_gamma_c = one_over_gamma*inv_c; + const auto bx = ux*one_over_gamma_c; + const auto by = uy*one_over_gamma_c; + const auto bz = uz*one_over_gamma_c; + + const auto one_over_dt_gamma_c = one_over_gamma_c/dt; + const auto bpx = (p_ux[ip] - p_ux_old[ip])*one_over_dt_gamma_c; + const auto bpy = (p_uy[ip] - p_uy_old[ip])*one_over_dt_gamma_c; + const auto bpz = (p_uz[ip] - p_uz_old[ip])*one_over_dt_gamma_c; + + const auto tot_q = q*p_w[ip]; + + for(int i_om=0; i_om < omega_points; ++i_om){ + + const auto i_omega_over_c = Complex{0.0_prt, 1.0_prt}*p_omegas[i_om]*inv_c; + for (int i_det = 0; i_det < how_many_det_pos; ++i_det){ + + const auto part_det_x = xp - p_det_pos_x[i_det]; + const auto part_det_y = yp - p_det_pos_y[i_det]; + const auto part_det_z = zp - p_det_pos_z[i_det]; + const auto d_part_det = std::sqrt( + part_det_x*part_det_x + part_det_y*part_det_y + part_det_z*part_det_z); + const auto one_over_d_part_det = 1.0_prt/d_part_det; + + const auto nx = part_det_x*one_over_d_part_det; + const auto ny = part_det_y*one_over_d_part_det; + const auto nz = part_det_z*one_over_d_part_det; + + //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal + const auto one_minus_b_dot_n = 1.0_prt - (bx*nx + by*ny + bz*nz); + + const auto n_minus_beta_x = nx - bx; + const auto n_minus_beta_y = ny - by; + const auto n_minus_beta_z = nz - bz; + + //Calculation of nxbeta + const auto n_minus_beta_cross_bp_x = n_minus_beta_y*bpz - n_minus_beta_z*bpy; + const auto n_minus_beta_cross_bp_y = n_minus_beta_z*bpx - n_minus_beta_x*bpz; + const auto n_minus_beta_cross_bp_z = n_minus_beta_x*bpy - n_minus_beta_y*bpx; + + //Calculation of nxnxbeta + const auto n_cross_n_minus_beta_cross_bp_x = ny*n_minus_beta_cross_bp_z - nz*n_minus_beta_cross_bp_y; + const auto n_cross_n_minus_beta_cross_bp_y = nz*n_minus_beta_cross_bp_x - nx*n_minus_beta_cross_bp_z; + const auto n_cross_n_minus_beta_cross_bp_z = nx*n_minus_beta_cross_bp_y - ny*n_minus_beta_cross_bp_x; + + const auto phase_term = amrex::exp(i_omega_over_c*(c*current_time - d_part_det)); + + const auto coeff = tot_q*phase_term/(one_minus_b_dot_n*one_minus_b_dot_n); + + auto cx = coeff*n_cross_n_minus_beta_cross_bp_x; + auto cy = coeff*n_cross_n_minus_beta_cross_bp_y; + auto cz = coeff*n_cross_n_minus_beta_cross_bp_z; + + // Nyquist limiter + /*if(p_omegas[i_om] < ablastr::constant::math::pi/one_minus_b_dot_n/dt){ + cx = 0.0; + cy = 0.0; + cz = 0.0; + }*/ + + const int ncomp = 3; + const int idx0 = (i_om*how_many_det_pos + i_det)*ncomp; + const int idx1 = idx0 + 1; + const int idx2 = idx0 + 2; + +#ifdef AMREX_USE_OMP + #pragma omp atomic + p_radiation_data[idx0].m_real += cx.m_real; + #pragma omp atomic + p_radiation_data[idx0].m_imag += cx.m_imag; + #pragma omp atomic + p_radiation_data[idx1].m_real += cy.m_real; + #pragma omp atomic + p_radiation_data[idx1].m_imag += cy.m_imag; + #pragma omp atomic + p_radiation_data[idx2].m_real += cz.m_real; + #pragma omp atomic + p_radiation_data[idx2].m_imag += cz.m_imag; +#else + p_radiation_data[idx0] += cx; + p_radiation_data[idx1] += cy; + p_radiation_data[idx2] += cz; +#endif + } + } + }); + + } + } + } + } + + +void RadiationHandler::gather_and_write_radiation(const std::string& filename) +{ + auto radiation_data_cpu = amrex::Vector(m_det_pts[0]*m_det_pts[1]*m_omega_points); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + m_radiation_calculation.begin(), m_radiation_calculation.end(), radiation_data_cpu.begin()); + amrex::Gpu::streamSynchronize(); + + amrex::ParallelDescriptor::ReduceRealSum(radiation_data_cpu.data(), radiation_data_cpu.size()); + + if (amrex::ParallelDescriptor::IOProcessor() ){ + auto of = std::ofstream(filename, std::ios::binary); + + const auto how_many = m_det_pts[0]*m_det_pts[1]; + auto det_pos_x_cpu = amrex::Vector(how_many); + auto det_pos_y_cpu = amrex::Vector(how_many); + auto det_pos_z_cpu = amrex::Vector(how_many); + auto det_pos_theta_cpu = amrex::Vector(how_many); + auto det_pos_phi_cpu = amrex::Vector(how_many); + auto omegas_cpu = amrex::Vector(m_omega_points); + + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + det_pos_x.begin(), det_pos_x.end(), det_pos_x_cpu.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + det_pos_y.begin(), det_pos_y.end(), det_pos_y_cpu.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + det_pos_z.begin(), det_pos_z.end(), det_pos_z_cpu.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + det_pos_theta.begin(), det_pos_theta.end(), det_pos_theta_cpu.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + det_pos_phi.begin(), det_pos_phi.end(), det_pos_phi_cpu.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + m_omegas.begin(), m_omegas.end(), omegas_cpu.begin()); + amrex::Gpu::streamSynchronize(); + + + int idx = 0; + for(int i_om=0; i_om < m_omega_points; ++i_om){ + for (int i_det = 0; i_det < how_many; ++i_det) + { + of << omegas_cpu[i_om] << " " << det_pos_theta_cpu[i_det] << " " << det_pos_phi_cpu[i_det] << " " << det_pos_x_cpu[i_det] << " " << det_pos_y_cpu[i_det] << " " << det_pos_z_cpu[i_det] << " " << radiation_data_cpu[++idx] << "\n"; + } + } + + of.close(); + } +} + +void RadiationHandler::Integral_overtime(const amrex::Real dt) +{ + const auto factor = ablastr::constant::SI::q_e*dt/16/std::pow(ablastr::constant::math::pi,3)/PhysConst::ep0/(PhysConst::c); + const auto how_many = m_det_pts[0]*m_det_pts[1]; + auto p_radiation_data = m_radiation_data.dataPtr(); + m_radiation_calculation.resize(how_many*m_omega_points); + for(int idx=0; idx& v ); + /** \brief Find the indices that would reorder the elements of `predicate` * so that the elements with non-zero value precede the other elements * diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 024cd8d9157..71c202a342f 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -375,6 +375,8 @@ public: int DoQED() const { return false; } #endif + virtual bool has_radiation() const {return false;} + /* \brief This function tests if the current species * is of a given PhysicalSpecies (specified as a template parameter). * @tparam PhysSpec the PhysicalSpecies to test against @@ -404,7 +406,21 @@ public: */ void defineAllParticleTiles () noexcept; + /** + * This function copies the momenta of the particles into: + * old_u_x, old_u_y, old_u_z. + */ + void RecordOldMomenta(); + protected: + + /** + * This function allocates three real components: + * old_u_x, old_u_y, old_u_z. + */ + void AllocateOldMomenta(); + + int species_id; amrex::ParticleReal charge; @@ -430,6 +446,7 @@ protected: int do_continuous_injection = 0; int do_field_ionization = 0; + int ionization_product; std::string ionization_product_name; int ion_atomic_number; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 447721a3e68..76a74da325f 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1067,7 +1067,6 @@ WarpXParticleContainer::DepositCharge (std::unique_ptr& rho, } #endif } - std::unique_ptr WarpXParticleContainer::GetChargeDensity (int lev, bool local) { @@ -1416,3 +1415,44 @@ WarpXParticleContainer::ApplyBoundaryConditions (){ } } } + +void WarpXParticleContainer::AllocateOldMomenta() +{ + AddRealComp("old_u_x"); + AddRealComp("old_u_y"); + AddRealComp("old_u_z"); +} + +void WarpXParticleContainer::RecordOldMomenta() +{ + auto const nlevs = std::max(0, finestLevel()+1); + for (int lev = 0; lev < nlevs; ++lev) { + +#ifdef AMREX_USE_OMP + #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) +#endif + { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + auto& part= GetParticles(lev)[index]; + auto const np = pti.numParticles(); + auto& soa = part.GetStructOfArrays(); + + const auto* ux = soa.GetRealData(PIdx::ux).data(); + const auto* uy = soa.GetRealData(PIdx::uy).data(); + const auto* uz = soa.GetRealData(PIdx::uz).data(); + + auto* old_ux = soa.GetRealData(GetRealCompIndex("old_u_x")).data(); + auto* old_uy = soa.GetRealData(GetRealCompIndex("old_u_y")).data(); + auto* old_uz = soa.GetRealData(GetRealCompIndex("old_u_z")).data(); + + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE (int ip) { + old_ux[ip] = ux[ip]; + old_uy[ip] = uy[ip]; + old_uz[ip] = uz[ip]; + }); + + } + } + } +} diff --git a/Source/Utils/WarpX_Complex.H b/Source/Utils/WarpX_Complex.H deleted file mode 100644 index 7a1f1669286..00000000000 --- a/Source/Utils/WarpX_Complex.H +++ /dev/null @@ -1,32 +0,0 @@ -/* Copyright 2019-2020 Andrew Myers, David Grote, Maxence Thevenet - * Remi Lehe - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef WARPX_COMPLEX_H_ -#define WARPX_COMPLEX_H_ - -#ifdef WARPX_USE_PSATD -# include -#endif - -#include -#include -#include - -#include - -// Defines a complex type on GPU & CPU -using Complex = amrex::GpuComplex; - -#ifdef WARPX_USE_PSATD -static_assert(sizeof(Complex) == sizeof(ablastr::math::anyfft::Complex), - "The complex type in WarpX and the FFT library do not match."); -#endif - -static_assert(sizeof(Complex) == sizeof(amrex::Real[2]), - "Unexpected complex type."); - -#endif //WARPX_COMPLEX_H_ diff --git a/Source/ablastr/math/Complex.H b/Source/ablastr/math/Complex.H new file mode 100644 index 00000000000..705e5182bbb --- /dev/null +++ b/Source/ablastr/math/Complex.H @@ -0,0 +1,36 @@ +/* Copyright 2019-2020 Andrew Myers, David Grote, Maxence Thevenet + * Remi Lehe + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_MATH_COMPLEX_H_ +#define ABLASTR_MATH_COMPLEX_H_ + +#ifdef ABLASTR_USE_FFT +# include "ablastr/math/fft/AnyFFT.H" +#endif +#include "ablastr/utils/TextMsg.H" + +#include +#include +#include + +#include + +namespace ablastr::math +{ + // Defines a complex type on GPU & CPU + using Complex = amrex::GpuComplex; + +#ifdef ABLASTR_USE_FFT + static_assert(sizeof(Complex) == sizeof(ablastr::math::anyfft::Complex), + "The complex type in ablastr and the FFT library do not match."); +#endif + + static_assert(sizeof(Complex) == sizeof(amrex::Real[2]), + "Unexpected complex type."); +} + +#endif //ABLASTR_MATH_COMPLEX_H_ diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index 2b044f3a485..8d7e468b029 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -269,7 +269,7 @@ set(WarpX_amrex_src "" set(WarpX_amrex_repo "https://github.com/AMReX-Codes/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(WarpX_amrex_internal)") -set(WarpX_amrex_branch "ecaa46d0be4b5c79b8806e48e3469000d8bb7252" +set(WarpX_amrex_branch "f692e78ac94a39bb5a5e40c5019dac08d702f6cc" CACHE STRING "Repository branch for WarpX_amrex_repo if(WarpX_amrex_internal)")