From 9df2a10d47e1b8b42a8071d0ada4f580e5478097 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Mon, 4 Dec 2023 14:53:38 +0900 Subject: [PATCH 01/60] first_step --- Source/Evolve/WarpXEvolve.cpp | 12 ++ Source/Particles/CMakeLists.txt | 1 + Source/Particles/Make.package | 1 + Source/Particles/MultiParticleContainer.H | 7 +- Source/Particles/MultiParticleContainer.cpp | 24 +++- Source/Particles/PhysicalParticleContainer.H | 3 + .../Particles/PhysicalParticleContainer.cpp | 10 ++ Source/Particles/Radiation/CMakeLists.txt | 7 + Source/Particles/Radiation/Make.package | 1 + Source/Particles/Radiation/RadiationHandler.H | 66 +++++++++ .../Particles/Radiation/RadiationHandler.cpp | 136 ++++++++++++++++++ Source/Particles/WarpXParticleContainer.H | 3 + Source/Utils/Math/Complex.H | 56 ++++++++ Source/WarpX.H | 4 + Source/WarpX.cpp | 4 +- 15 files changed, 331 insertions(+), 4 deletions(-) create mode 100644 Source/Particles/Radiation/CMakeLists.txt create mode 100644 Source/Particles/Radiation/Make.package create mode 100644 Source/Particles/Radiation/RadiationHandler.H create mode 100644 Source/Particles/Radiation/RadiationHandler.cpp create mode 100644 Source/Utils/Math/Complex.H diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 1515a1b3ff0..63be9cd8d0a 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -427,6 +427,18 @@ WarpX::OneStep_nosub (Real cur_time) PushParticlesandDepose(cur_time); + //Radiation contribution at each timestep + if (m_do_radiation_flag){ + //Save particle old momentum in a attribute + //mypc->RadiationHandler::keepoldmomentum(); + //Only level 0 is supported + mypc->doRadiation(dt[0]); + } + + + + + ExecutePythonCallback("afterdeposition"); // Synchronize J and rho: 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/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 64ce057e2ff..10e614b8c0a 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.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" #include "WarpXParticleContainer.H" @@ -329,6 +330,8 @@ public: const amrex::MultiFab& Bz); #endif + void doRadiation(amrex::Real dt); + int getSpeciesID (std::string product_str) const; protected: @@ -526,6 +529,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 5cfd4e2b68d..064798251ed 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -388,8 +388,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); @@ -404,6 +403,16 @@ MultiParticleContainer::ReadParameters () #endif initialized = true; } + + + for (auto& s : allcontainers) { + + if (s->has_radiation()) { + m_at_least_one_has_radiation = true; + m_p_radiation_handler = std::make_unique(); + break; + } + } } WarpXParticleContainer& @@ -975,6 +984,17 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vectorhas_radiation()){ + m_p_radiation_handler->add_radiation_contribution(dt,pc); + + } + } +} + #ifdef WARPX_QED void MultiParticleContainer::InitQED () { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 526e5e01b2d..cb60490700d 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -353,6 +353,8 @@ public: PairGenerationFilterFunc getPairGenerationFilterFunc (); #endif + bool has_radiation() const override {return m_do_radiation;} + protected: std::string species_name; std::unique_ptr plasma_injector; @@ -381,6 +383,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 47898a0344d..3f059b4f448 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -279,6 +279,16 @@ 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 previous momentum attribute + pp_species_name.query("do_radiation", m_do_radiation); + if (m_do_radiation){ + //enable the radiations in the evolve looop + WarpX::m_do_radiation_flag = 1; + + AddRealComp("previous_momentum_x"); + AddRealComp("previous_momentum_y"); + AddRealComp("previous_momentum_z"); + } //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..70659a0e305 --- /dev/null +++ b/Source/Particles/Radiation/Make.package @@ -0,0 +1 @@ +CEXE_sources += RadiationHandler.cpp diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H new file mode 100644 index 00000000000..b0e49cf5d6f --- /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 "Utils/Math/Complex.H" +#include "Particles/WarpXParticleContainer.H" +#include "Particles/Pusher/GetAndSetPosition.H" + + +#include +#include +#include +#include +#include + +#include +#include + +/* \brief CollisionHandler is a light weight class that contains the + * calculation of radiation for particles at each frequencies and angles. + */ +class RadiationHandler +{ +public: + RadiationHandler (); + + /* 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); + void add_detector(amrex::IntVect m_det_pts, amrex::IntVect m_omega_rg, amrex::Vector m_d_det, int m_omega_pts, double m_d_distance, amrex::IntVect m_d_direction, amrex::Vector m_th_rg); + void keepoldmomentum(std::unique_ptr& pc); + +private: + + // Frequency range + amrex::IntVect m_omega_range; + + // Dimensions of the detector + amrex::Vector m_theta_range; + + int m_omega_points; + + //put a detector + int m_get_a_detector; + // Resolution of the detector + amrex::Vector m_d_theta; + double m_d_omega; + //The vector with the resolution of the detector + amrex::Vector m_d_det; + int m_det_distance; + + //Define the Fab with the datas + int ncomp; + + amrex::IntVect m_det_pts; + amrex::IntVect m_det_direction; + amrex::Vector m_pos_det; +}; +#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..3bdc7f05d8b --- /dev/null +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -0,0 +1,136 @@ +/* Copyright 2023 Thomas Clark, Luca Fedeli + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "RadiationHandler.H" + +#include "Utils/TextMsg.H" + +#include + +#include +using namespace warpx::utils; + +RadiationHandler::RadiationHandler() +{ + // Read in radiation input + const amrex::ParmParse pp_radiations("radiations"); + + // Direction of the normal of the detector + pp_radiations.query("omega_range", m_omega_range); + pp_radiations.query("omega_points", m_omega_points); + + pp_radiations.queryarr("angle_aperture", m_theta_range); + + pp_radiations.query("put_a_detector", m_get_a_detector); + + +#if defined WARPX_DIM_RZ + WARPX_ABORT_WITH_MESSAGE("Radiation is not supported yet with RZ."); +#endif + + //Compute the radiation + if(m_get_a_detector==1){ + pp_radiations.query("detector_number_points", m_det_pts); + pp_radiations.query("detector_direction", m_det_direction); + pp_radiations.query("detector_distance", m_det_distance); + + add_detector(m_det_pts, m_omega_range, m_d_det, m_omega_points, m_det_distance, m_det_direction, m_theta_range); + } + + // + + + /* Initialize the Fab with the field in function of angle and frequency */ + //int numcomps = 4; + //BaseFab fab(bx,numcomps); +} + +void RadiationHandler::add_radiation_contribution + (const amrex::Real dt, std::unique_ptr& pc){ + const auto level0=0; + const auto part=pc->GetParticles(level0); + +#ifdef AMREX_USE_OMP +#pragma omp parallel for +#endif + for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { + long const np = pti.numParticles(); + auto& attribs = pti.GetAttribs(); + auto& p_w = attribs[PIdx::w]; + auto& p_ux = attribs[PIdx::ux]; + auto& p_uy = attribs[PIdx::uy]; + auto& p_uz = attribs[PIdx::uz]; + auto& attribut = attribs[PIdx::nattribs]; + + //auto GetPosition = GetParticlePosition(pti); + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int ip) + { //amrex::ParticleReal xp, yp, zp; + //GetPosition.AsStored(ip, xp, yp, zp); + const amrex::ParticleReal p_ux_unit=p_ux[ip]; + const amrex::ParticleReal p_uy_unit=p_uy[ip]; + const amrex::ParticleReal p_uz_unit=p_uz[ip]; + + amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); + + //const amrex::ParticleReal* p_pos0 = src_data.pos(0); + + }); + +// const int total_partdiag_size = amrex::Scan::ExclusiveSum(np,Flag,IndexLocation); +// auto& ptile_dst = pc_dst.DefineAndReturnParticleTile(lev, pti.index(), pti.LocalTileIndex() ); +// auto old_size = ptile_dst.numParticles(); +// ptile_dst.resize(old_size + total_partdiag_size); +// amrex::filterParticles(ptile_dst, ptile_src, GetParticleFilter, 0, old_size, np); +// auto dst_data = ptile_dst.getParticleTileData(); +// amrex::ParallelFor(np, +// [=] AMREX_GPU_DEVICE(int i) +// { +// if (Flag[i] == 1) GetParticleLorentzTransform(dst_data, src_data, i, +// old_size + IndexLocation[i]); +// }); +// amrex::Gpu::synchronize(); + } + +} + +void RadiationHandler::add_detector + (amrex::IntVect m_d_pts, amrex::IntVect m_omega_rg, amrex::Vector m_d_d, int m_omega_pts, double m_d_distance, amrex::IntVect m_d_direction, amrex::Vector m_th_rg){ + + // Box of angle and frequency + const amrex::Box detect_box({0,0,0}, {m_d_pts[0], m_d_pts[1], m_omega_pts}); + amrex::Print() << detect_box << "\n"; + ncomp = 0; + amrex::FArrayBox fab_detect(detect_box, ncomp); + amrex::Print() << fab_detect << "\n"; + + //Calculation of angle resolution + for(int i=0; i<2; i++){ + m_d_theta[i] = 2*m_th_rg[i]/static_cast(m_d_pts[i]); + } + //Set the resolution of the detector + for(int idim = 0; idim<3; ++idim){ + if(m_d_direction[idim]==1){ + m_d_d[(idim+1)%3]=2*m_d_distance*tan(m_d_theta[0]/2); + m_d_d[(idim+2)%3]=2*m_d_distance*tan(m_d_theta[1]/2); + } + + m_d_omega=m_omega_rg[1]-m_omega_rg[0]/static_cast(m_omega_pts); + } +} + +void RadiationHandler::keepoldmomentum + (std::unique_ptr& pc){ + const auto level0=0; + const auto part=pc->GetParticles(level0); + + for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { + //long const np = pti.numParticles(); + //auto& attribs = pti.GetAtt + }; + +} diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 8587f74544d..0e86ff6931c 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -366,6 +366,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 @@ -421,6 +423,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/Utils/Math/Complex.H b/Source/Utils/Math/Complex.H new file mode 100644 index 00000000000..81a29407a31 --- /dev/null +++ b/Source/Utils/Math/Complex.H @@ -0,0 +1,56 @@ +/* Copyright 2022 Thomas Clark + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_UTILS_MATH_COMPLEX__H_ +#define WARPX_UTILS_MATH_COMPLEX__H_ + +#include +#include + +#if defined(AMREX_USE_CUDA) +# include +#elif defined(AMREX_USE_HIP) +# include + +#else +# if WARPX_USE_PSATD +# include +# else +# include + +# endif +#endif + + // First, define library-dependent types (complex, FFT plan) +namespace warpx::utils +{ + /** Complex type for FFT, depends on FFT library */ +#if defined(AMREX_USE_CUDA) +# ifdef AMREX_USE_FLOAT + using Complex = cuComplex; +# else + using Complex = cuDoubleComplex; +# endif +#elif defined(AMREX_USE_HIP) +# ifdef AMREX_USE_FLOAT + using Complex = float2; +# else + using Complex = double2; +# endif +#else +#if WARPX_USE_PSATD +# ifdef AMREX_USE_FLOAT + using Complex = fftwf_complex; +# else + using Complex = fftw_complex; +# endif +# else + using Complex = std::complex; +# endif +#endif +} +#endif // WARPX_UTILS_MATH_COMPLEX__H_ \ No newline at end of file diff --git a/Source/WarpX.H b/Source/WarpX.H index 6a4d76f3c01..aa3630e183a 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -354,6 +354,10 @@ public: static bool do_device_synchronize; static bool safe_guard_cells; + + //enable radiations in the evolve + static bool m_do_radiation_flag; + //! With mesh refinement, particles located inside a refinement patch, but within //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields //! from the lower refinement level instead of the refinement patch itself diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 240c7fbb4e0..ad135820d3e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -222,6 +222,8 @@ bool WarpX::do_device_synchronize = true; bool WarpX::do_device_synchronize = false; #endif +bool m_do_radiation_flag = 0; + WarpX* WarpX::m_instance = nullptr; void WarpX::MakeWarpX () @@ -1928,7 +1930,7 @@ WarpX::BackwardCompatibility () "collisions.ncollisions is ignored. Just use particles.collision_names please.", ablastr::warn_manager::WarnPriority::low); } - + const ParmParse pp_lasers("lasers"); int nlasers; if (pp_lasers.query("nlasers", nlasers)){ From 3fc03cb0e134273685ff246923b025ec81b11697 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Mon, 4 Dec 2023 18:50:07 +0900 Subject: [PATCH 02/60] update --- .../Particles/PhysicalParticleContainer.cpp | 6 +++--- .../Particles/Radiation/RadiationHandler.cpp | 19 ++++++++++++------- 2 files changed, 15 insertions(+), 10 deletions(-) diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3f059b4f448..63fc7fe31cc 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -285,9 +285,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //enable the radiations in the evolve looop WarpX::m_do_radiation_flag = 1; - AddRealComp("previous_momentum_x"); - AddRealComp("previous_momentum_y"); - AddRealComp("previous_momentum_z"); + AddRealComp("prev_u_x"); + AddRealComp("prev_u_y"); + AddRealComp("prev_u_z"); } //if the species is not a lepton, do_classical_radiation_reaction //should be false diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 3bdc7f05d8b..dbfe862ed7b 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -57,6 +57,7 @@ void RadiationHandler::add_radiation_contribution #ifdef AMREX_USE_OMP #pragma omp parallel for #endif +{ for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { long const np = pti.numParticles(); auto& attribs = pti.GetAttribs(); @@ -66,6 +67,10 @@ void RadiationHandler::add_radiation_contribution auto& p_uz = attribs[PIdx::uz]; auto& attribut = attribs[PIdx::nattribs]; + //auto& ux = src.m_rdata[PIdx::ux][i_src]; + //auto& uy = src.m_rdata[PIdx::uy][i_src]; + //auto& uz = src.m_rdata[PIdx::uz][i_src]; + //auto GetPosition = GetParticlePosition(pti); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip) @@ -94,8 +99,8 @@ void RadiationHandler::add_radiation_contribution // old_size + IndexLocation[i]); // }); // amrex::Gpu::synchronize(); + } } - } void RadiationHandler::add_detector @@ -126,11 +131,11 @@ void RadiationHandler::add_detector void RadiationHandler::keepoldmomentum (std::unique_ptr& pc){ const auto level0=0; - const auto part=pc->GetParticles(level0); - for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { - //long const np = pti.numParticles(); - //auto& attribs = pti.GetAtt - }; - + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + auto& part=pc->GetParticles(level0)[index]; + long const np = pti.numParticles(); + auto& soa = part.GetStructOfArrays(); + const amrex::ParticleReal* const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + pc->particle_runtime_comps["prev_u_x"]; } From a9728a584cd6e3587796e236a239201acb7734d2 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Tue, 5 Dec 2023 10:12:38 +0900 Subject: [PATCH 03/60] . --- Source/Evolve/WarpXEvolve.cpp | 9 +++++---- Source/Particles/MultiParticleContainer.H | 8 ++++++++ Source/Particles/MultiParticleContainer.cpp | 15 ++++++++++++++- .../Particles/NamedComponentParticleContainer.H | 17 +++++++++++++++++ Source/Particles/PhysicalParticleContainer.cpp | 2 -- Source/Particles/Radiation/RadiationHandler.H | 1 - Source/Particles/Radiation/RadiationHandler.cpp | 13 ++++++++++--- Source/WarpX.H | 2 +- 8 files changed, 55 insertions(+), 12 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 63be9cd8d0a..f05cb31e48c 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -424,16 +424,17 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("particlescraper"); ExecutePythonCallback("beforedeposition"); + if (m_do_radiation_flag){ + //Save particle old momentum in a attribute + //mypc->RadiationHandler::keepoldmomentum(); + } PushParticlesandDepose(cur_time); //Radiation contribution at each timestep - if (m_do_radiation_flag){ - //Save particle old momentum in a attribute - //mypc->RadiationHandler::keepoldmomentum(); //Only level 0 is supported mypc->doRadiation(dt[0]); - } + diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 10e614b8c0a..dcdadb4717c 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -332,6 +332,8 @@ public: void doRadiation(amrex::Real dt); + void initradiation(); + int getSpeciesID (std::string product_str) const; protected: @@ -482,6 +484,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: diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 064798251ed..860ea004c72 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -983,14 +983,27 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vectorhas_radiation()){ + m_at_least_one_particle_radiate = 1; + m_nspecies_radiate++; + } + } + +} void MultiParticleContainer::doRadiation (const amrex::Real dt) -{ +{ if (m_at_least_one_has_radiation){ for (auto& pc : allcontainers) { if (pc->has_radiation()){ m_p_radiation_handler->add_radiation_contribution(dt,pc); + } } } } diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index 6ff9d157790..7986cb0cc31 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -176,6 +176,23 @@ public: std::map getParticleRuntimeComps () const noexcept { return particle_runtime_comps;} /** Return the name-to-index map for the runtime-time integer components */ 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_runtime_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; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 63fc7fe31cc..48985770415 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -283,8 +283,6 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp_species_name.query("do_radiation", m_do_radiation); if (m_do_radiation){ //enable the radiations in the evolve looop - WarpX::m_do_radiation_flag = 1; - AddRealComp("prev_u_x"); AddRealComp("prev_u_y"); AddRealComp("prev_u_z"); diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index b0e49cf5d6f..2dd62365899 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -9,7 +9,6 @@ #ifndef WARPX_PARTICLES_RADIATION_H_ #define WARPX_PARTICLES_RADIATION_H_ #include "Utils/Math/Complex.H" -#include "Particles/WarpXParticleContainer.H" #include "Particles/Pusher/GetAndSetPosition.H" diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index dbfe862ed7b..07b10a827ba 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -136,6 +136,13 @@ void RadiationHandler::keepoldmomentum auto& part=pc->GetParticles(level0)[index]; long const np = pti.numParticles(); auto& soa = part.GetStructOfArrays(); - const amrex::ParticleReal* const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); - pc->particle_runtime_comps["prev_u_x"]; -} + amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); + int index_name=pc->GetRealCompIndex("prev_u_x"); + auto& p_ux = soa.GetRealData(index_name); + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int ip) + { + p_ux[ip] = ux[ip]; + }); + } + } diff --git a/Source/WarpX.H b/Source/WarpX.H index aa3630e183a..e28d859ffd4 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -356,7 +356,7 @@ public: //enable radiations in the evolve - static bool m_do_radiation_flag; + bool m_do_radiation_flag; //! With mesh refinement, particles located inside a refinement patch, but within //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields From 2f3384e76b160f7fac35dd0f8347e77ed3b949c1 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Tue, 5 Dec 2023 11:36:21 +0900 Subject: [PATCH 04/60] . --- Source/Evolve/WarpXEvolve.cpp | 5 ++- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 28 ++++++--------- .../Particles/PhysicalParticleContainer.cpp | 2 +- Source/Particles/Radiation/RadiationHandler.H | 1 - .../Particles/Radiation/RadiationHandler.cpp | 21 +----------- Source/Particles/WarpXParticleContainer.H | 3 ++ Source/Particles/WarpXParticleContainer.cpp | 34 +++++++++++++++++++ 8 files changed, 53 insertions(+), 43 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index f05cb31e48c..e7e910cef59 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -424,10 +424,9 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("particlescraper"); ExecutePythonCallback("beforedeposition"); - if (m_do_radiation_flag){ + //Save particle old momentum in a attribute - //mypc->RadiationHandler::keepoldmomentum(); - } + mypc->keepoldmomentum(); PushParticlesandDepose(cur_time); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index dcdadb4717c..6aba6d20d81 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -332,7 +332,7 @@ public: void doRadiation(amrex::Real dt); - void initradiation(); + void keepoldmomentum(); int getSpeciesID (std::string product_str) const; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 860ea004c72..371bd41d874 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -404,7 +404,7 @@ MultiParticleContainer::ReadParameters () initialized = true; } - +//Initialization of the radiation for (auto& s : allcontainers) { if (s->has_radiation()) { @@ -983,32 +983,26 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vectorhas_radiation()){ - m_at_least_one_particle_radiate = 1; - m_nspecies_radiate++; + if (pc->has_radiation()){ + //m_p_radiation_handler->add_radiation_contribution(dt,pc); + + } } } - } - -void -MultiParticleContainer::doRadiation (const amrex::Real dt) -{ if (m_at_least_one_has_radiation){ +#ifdef WARPX_QED +void MultiParticleContainer::keepoldmomentum(){ for (auto& pc : allcontainers) { if (pc->has_radiation()){ - m_p_radiation_handler->add_radiation_contribution(dt,pc); - - } + pc->keepoldmomentum(); } } } - -#ifdef WARPX_QED void MultiParticleContainer::InitQED () { m_shr_p_qs_engine = std::make_shared(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 48985770415..be91d00c5e6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -282,7 +282,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //check if Radiations are enables for species and add previous momentum attribute pp_species_name.query("do_radiation", m_do_radiation); if (m_do_radiation){ - //enable the radiations in the evolve looop + //enable the radiations in the evolve loop AddRealComp("prev_u_x"); AddRealComp("prev_u_y"); AddRealComp("prev_u_z"); diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 2dd62365899..61982d155de 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -33,7 +33,6 @@ public: //void doRadiation (amrex::Real dt, MultiParticleContainer* mypc); void add_radiation_contribution(const amrex::Real dt, std::unique_ptr& pc); void add_detector(amrex::IntVect m_det_pts, amrex::IntVect m_omega_rg, amrex::Vector m_d_det, int m_omega_pts, double m_d_distance, amrex::IntVect m_d_direction, amrex::Vector m_th_rg); - void keepoldmomentum(std::unique_ptr& pc); private: diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 07b10a827ba..9006a9412c0 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -126,23 +126,4 @@ void RadiationHandler::add_detector m_d_omega=m_omega_rg[1]-m_omega_rg[0]/static_cast(m_omega_pts); } -} - -void RadiationHandler::keepoldmomentum - (std::unique_ptr& pc){ - const auto level0=0; - for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { - auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); - auto& part=pc->GetParticles(level0)[index]; - long const np = pti.numParticles(); - auto& soa = part.GetStructOfArrays(); - amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); - int index_name=pc->GetRealCompIndex("prev_u_x"); - auto& p_ux = soa.GetRealData(index_name); - amrex::ParallelFor(np, - [=] AMREX_GPU_DEVICE(int ip) - { - p_ux[ip] = ux[ip]; - }); - } - } +} \ No newline at end of file diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 0e86ff6931c..f257aefa11f 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -217,6 +217,9 @@ public: bool apply_boundary_and_scale_volume = false, int icomp = 0); + //Keep old momentum for particles + void keepoldmomentum (std::unique_ptr& pc); + std::unique_ptr GetChargeDensity(int lev, bool local = false); virtual void DepositCharge (WarpXParIter& pti, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 6109cd830f9..ea9a6efaef8 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1055,6 +1055,40 @@ WarpXParticleContainer::DepositCharge (std::unique_ptr& rho, } #endif } +void WarpXParticleContainer::keepoldmomentum + (std::unique_ptr& pc){ + const auto level0=0; + for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + auto& part=pc->GetParticles(level0)[index]; + long const np = pti.numParticles(); + auto& soa = part.GetStructOfArrays(); + + //Load the momentums + amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); + amrex::ParticleReal* uy = soa.GetRealData(PIdx::ux).data(); + amrex::ParticleReal* uz = soa.GetRealData(PIdx::ux).data(); + + //Finding the good attribute index + int index_name_x=pc->GetRealCompIndex("prev_u_x"); + int index_name_y=pc->GetRealCompIndex("prev_u_y"); + int index_name_z=pc->GetRealCompIndex("prev_u_z"); + + auto* p_ux = soa.GetRealData(index_name_x).data(); + auto* p_uy = soa.GetRealData(index_name_y).data(); + auto* p_uz = soa.GetRealData(index_name_z).data(); + + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int ip) + { + //Putting them in an attribute + p_ux[ip] = ux[ip]; + p_uy[ip] = ux[ip]; + p_uz[ip] = ux[ip]; + + }); + } + } std::unique_ptr WarpXParticleContainer::GetChargeDensity (int lev, bool local) From f52b7f6cc0d64b5936d93a65d16023b5081a577e Mon Sep 17 00:00:00 2001 From: tmsclark Date: Tue, 5 Dec 2023 14:31:03 +0900 Subject: [PATCH 05/60] . --- Source/Particles/MultiParticleContainer.H | 2 + Source/Particles/MultiParticleContainer.cpp | 55 +++++++++++++++---- .../Particles/PhysicalParticleContainer.cpp | 3 +- Source/Particles/Radiation/RadiationHandler.H | 6 +- .../Particles/Radiation/RadiationHandler.cpp | 32 ++++++----- Source/Particles/WarpXParticleContainer.H | 3 - Source/Particles/WarpXParticleContainer.cpp | 35 ------------ 7 files changed, 71 insertions(+), 65 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 6aba6d20d81..a931a00f636 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -333,6 +333,8 @@ public: void doRadiation(amrex::Real dt); void keepoldmomentum(); + void keepoldmomentum_p(std::unique_ptr& pc); + int getSpeciesID (std::string product_str) const; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 371bd41d874..52912b04260 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -403,18 +403,18 @@ MultiParticleContainer::ReadParameters () #endif initialized = true; } - +amrex::Print() << "wEEEEEEEESSSSSSSHHHHH" << allcontainers.size() << "\n"; //Initialization of the radiation for (auto& s : allcontainers) { - + amrex::Print() << s->has_radiation() <<"jzfnfhjzbefhezbfhzebfzehb" << "\n"; if (s->has_radiation()) { m_at_least_one_has_radiation = true; + amrex::Print() << "Hey ! I'm here" << "\n"; m_p_radiation_handler = std::make_unique(); break; } } } - WarpXParticleContainer& MultiParticleContainer::GetParticleContainerFromName (const std::string& name) const { @@ -944,6 +944,47 @@ MultiParticleContainer::doFieldIonization (int lev, } } +void MultiParticleContainer::keepoldmomentum(){ + for (auto& pc : allcontainers) { + if (pc->has_radiation()){ + keepoldmomentum_p(pc); + } + } +}void MultiParticleContainer::keepoldmomentum_p + (std::unique_ptr& pc){ + const auto level0=0; + for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + auto& part=pc->GetParticles(level0)[index]; + long const np = pti.numParticles(); + auto& soa = part.GetStructOfArrays(); + + //Load the momentums + amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); + amrex::ParticleReal* uy = soa.GetRealData(PIdx::ux).data(); + amrex::ParticleReal* uz = soa.GetRealData(PIdx::ux).data(); + + //Finding the good attribute index + int index_name_x=pc->GetRealCompIndex("prev_u_x"); + int index_name_y=pc->GetRealCompIndex("prev_u_y"); + int index_name_z=pc->GetRealCompIndex("prev_u_z"); + + auto* p_ux = soa.GetRealData(index_name_x).data(); + auto* p_uy = soa.GetRealData(index_name_y).data(); + auto* p_uz = soa.GetRealData(index_name_z).data(); + + amrex::ParallelFor(np, + [=] AMREX_GPU_DEVICE(int ip) + { + //Putting them in an attribute + p_ux[ip] = ux[ip]; + p_uy[ip] = uy[ip]; + p_uz[ip] = uz[ip]; + + }); + } + } + void MultiParticleContainer::doCollisions ( Real cur_time, amrex::Real dt ) { @@ -989,6 +1030,7 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt) if (m_at_least_one_has_radiation){ for (auto& pc : allcontainers) { if (pc->has_radiation()){ + //m_p_radiation_handler->add_radiation_contribution(dt,pc); } @@ -996,13 +1038,6 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt) } } #ifdef WARPX_QED -void MultiParticleContainer::keepoldmomentum(){ - for (auto& pc : allcontainers) { - if (pc->has_radiation()){ - pc->keepoldmomentum(); - } - } -} void MultiParticleContainer::InitQED () { m_shr_p_qs_engine = std::make_shared(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index be91d00c5e6..b94c5d747b8 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -281,11 +281,12 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp_species_name.query("do_classical_radiation_reaction", do_classical_radiation_reaction); //check if Radiations are enables for species and add previous momentum attribute pp_species_name.query("do_radiation", m_do_radiation); - if (m_do_radiation){ + if (m_do_radiation==1){ //enable the radiations in the evolve loop AddRealComp("prev_u_x"); AddRealComp("prev_u_y"); AddRealComp("prev_u_z"); + } //if the species is not a lepton, do_classical_radiation_reaction //should be false diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 61982d155de..fe8ce97c576 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -32,7 +32,7 @@ public: /* 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); - void add_detector(amrex::IntVect m_det_pts, amrex::IntVect m_omega_rg, amrex::Vector m_d_det, int m_omega_pts, double m_d_distance, amrex::IntVect m_d_direction, amrex::Vector m_th_rg); + void add_detector(); private: @@ -59,6 +59,10 @@ private: amrex::IntVect m_det_pts; amrex::IntVect m_det_direction; amrex::Vector m_pos_det; + + //resolution of the grid of the detectir + amrex::Vector m_d_d; + }; #endif // WARPX_PARTICLES_RADIATION_H_ diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 9006a9412c0..889b8cc852f 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -19,12 +19,7 @@ RadiationHandler::RadiationHandler() // Read in radiation input const amrex::ParmParse pp_radiations("radiations"); - // Direction of the normal of the detector - pp_radiations.query("omega_range", m_omega_range); - pp_radiations.query("omega_points", m_omega_points); - - pp_radiations.queryarr("angle_aperture", m_theta_range); - + // Verify if there is a detector pp_radiations.query("put_a_detector", m_get_a_detector); @@ -33,12 +28,19 @@ RadiationHandler::RadiationHandler() #endif //Compute the radiation - if(m_get_a_detector==1){ + if(m_get_a_detector){ + //Resolution in frequency of the detector + pp_radiations.query("omega_range", m_omega_range); + pp_radiations.query("omega_points", m_omega_points); + //Angle theta AND phi + pp_radiations.queryarr("angle_aperture", m_theta_range); + + //Type of detector pp_radiations.query("detector_number_points", m_det_pts); pp_radiations.query("detector_direction", m_det_direction); pp_radiations.query("detector_distance", m_det_distance); - add_detector(m_det_pts, m_omega_range, m_d_det, m_omega_points, m_det_distance, m_det_direction, m_theta_range); + add_detector(); } // @@ -104,10 +106,10 @@ void RadiationHandler::add_radiation_contribution } void RadiationHandler::add_detector - (amrex::IntVect m_d_pts, amrex::IntVect m_omega_rg, amrex::Vector m_d_d, int m_omega_pts, double m_d_distance, amrex::IntVect m_d_direction, amrex::Vector m_th_rg){ + (){ // Box of angle and frequency - const amrex::Box detect_box({0,0,0}, {m_d_pts[0], m_d_pts[1], m_omega_pts}); + const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); amrex::Print() << detect_box << "\n"; ncomp = 0; amrex::FArrayBox fab_detect(detect_box, ncomp); @@ -115,15 +117,15 @@ void RadiationHandler::add_detector //Calculation of angle resolution for(int i=0; i<2; i++){ - m_d_theta[i] = 2*m_th_rg[i]/static_cast(m_d_pts[i]); + m_d_theta[i] = 2*m_theta_range[i]/static_cast(m_det_pts[i]); } //Set the resolution of the detector for(int idim = 0; idim<3; ++idim){ - if(m_d_direction[idim]==1){ - m_d_d[(idim+1)%3]=2*m_d_distance*tan(m_d_theta[0]/2); - m_d_d[(idim+2)%3]=2*m_d_distance*tan(m_d_theta[1]/2); + if(m_det_direction[idim]==1){ + m_d_d[(idim+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); + m_d_d[(idim+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); } - m_d_omega=m_omega_rg[1]-m_omega_rg[0]/static_cast(m_omega_pts); + m_d_omega=m_omega_range[1]-m_omega_range[0]/static_cast(m_omega_points); } } \ No newline at end of file diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index f257aefa11f..0e86ff6931c 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -217,9 +217,6 @@ public: bool apply_boundary_and_scale_volume = false, int icomp = 0); - //Keep old momentum for particles - void keepoldmomentum (std::unique_ptr& pc); - std::unique_ptr GetChargeDensity(int lev, bool local = false); virtual void DepositCharge (WarpXParIter& pti, diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index ea9a6efaef8..aaf6b4f48ab 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -1055,41 +1055,6 @@ WarpXParticleContainer::DepositCharge (std::unique_ptr& rho, } #endif } -void WarpXParticleContainer::keepoldmomentum - (std::unique_ptr& pc){ - const auto level0=0; - for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { - auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); - auto& part=pc->GetParticles(level0)[index]; - long const np = pti.numParticles(); - auto& soa = part.GetStructOfArrays(); - - //Load the momentums - amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); - amrex::ParticleReal* uy = soa.GetRealData(PIdx::ux).data(); - amrex::ParticleReal* uz = soa.GetRealData(PIdx::ux).data(); - - //Finding the good attribute index - int index_name_x=pc->GetRealCompIndex("prev_u_x"); - int index_name_y=pc->GetRealCompIndex("prev_u_y"); - int index_name_z=pc->GetRealCompIndex("prev_u_z"); - - auto* p_ux = soa.GetRealData(index_name_x).data(); - auto* p_uy = soa.GetRealData(index_name_y).data(); - auto* p_uz = soa.GetRealData(index_name_z).data(); - - amrex::ParallelFor(np, - [=] AMREX_GPU_DEVICE(int ip) - { - //Putting them in an attribute - p_ux[ip] = ux[ip]; - p_uy[ip] = ux[ip]; - p_uz[ip] = ux[ip]; - - }); - } - } - std::unique_ptr WarpXParticleContainer::GetChargeDensity (int lev, bool local) { From b6a3b1c53d3615042613da27c3d96aac722ad193 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Wed, 6 Dec 2023 13:53:35 +0900 Subject: [PATCH 06/60] . --- Source/Particles/MultiParticleContainer.cpp | 25 +++++++++---------- .../Particles/PhysicalParticleContainer.cpp | 2 +- Source/Particles/Radiation/RadiationHandler.H | 4 +-- .../Particles/Radiation/RadiationHandler.cpp | 17 ++++++------- 4 files changed, 23 insertions(+), 25 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 52912b04260..3bcf4a9781c 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -122,7 +122,17 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) // Setup particle collisions collisionhandler = std::make_unique(this); - + amrex::Print() << "wEEEEEEEESSSSSSSHHHHH" << allcontainers.size() << "\n"; + //Initialization of the radiation + for (auto& s : allcontainers) { + amrex::Print() << s->has_radiation() <<"jzfnfhjzbefhezbfhzebfzehb" << "\n"; + if (s->has_radiation()) { + m_at_least_one_has_radiation = true; + amrex::Print() << "Hey ! I'm here" << "\n"; + m_p_radiation_handler = std::make_unique(); + break; + } + } } void @@ -403,17 +413,6 @@ MultiParticleContainer::ReadParameters () #endif initialized = true; } -amrex::Print() << "wEEEEEEEESSSSSSSHHHHH" << allcontainers.size() << "\n"; -//Initialization of the radiation - for (auto& s : allcontainers) { - amrex::Print() << s->has_radiation() <<"jzfnfhjzbefhezbfhzebfzehb" << "\n"; - if (s->has_radiation()) { - m_at_least_one_has_radiation = true; - amrex::Print() << "Hey ! I'm here" << "\n"; - m_p_radiation_handler = std::make_unique(); - break; - } - } } WarpXParticleContainer& MultiParticleContainer::GetParticleContainerFromName (const std::string& name) const @@ -1031,7 +1030,7 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt) for (auto& pc : allcontainers) { if (pc->has_radiation()){ - //m_p_radiation_handler->add_radiation_contribution(dt,pc); + m_p_radiation_handler->add_radiation_contribution(dt,pc); } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b94c5d747b8..22500580d65 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -282,11 +282,11 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp //check if Radiations are enables for species and add previous momentum attribute pp_species_name.query("do_radiation", m_do_radiation); if (m_do_radiation==1){ + amrex::Print() << "i'm here"; //enable the radiations in the evolve loop AddRealComp("prev_u_x"); AddRealComp("prev_u_y"); AddRealComp("prev_u_z"); - } //if the species is not a lepton, do_classical_radiation_reaction //should be false diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index fe8ce97c576..ff05811625b 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -37,7 +37,7 @@ public: private: // Frequency range - amrex::IntVect m_omega_range; + amrex::Vector m_omega_range; // Dimensions of the detector amrex::Vector m_theta_range; @@ -56,7 +56,7 @@ private: //Define the Fab with the datas int ncomp; - amrex::IntVect m_det_pts; + amrex::Vector m_det_pts; amrex::IntVect m_det_direction; amrex::Vector m_pos_det; diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 889b8cc852f..ad511ceb674 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -30,16 +30,15 @@ RadiationHandler::RadiationHandler() //Compute the radiation if(m_get_a_detector){ //Resolution in frequency of the detector - pp_radiations.query("omega_range", m_omega_range); + pp_radiations.queryarr("omega_range", m_omega_range); pp_radiations.query("omega_points", m_omega_points); //Angle theta AND phi pp_radiations.queryarr("angle_aperture", m_theta_range); //Type of detector - pp_radiations.query("detector_number_points", m_det_pts); + pp_radiations.getarr("detector_number_points", m_det_pts,0,2); pp_radiations.query("detector_direction", m_det_direction); pp_radiations.query("detector_distance", m_det_distance); - add_detector(); } @@ -63,12 +62,12 @@ void RadiationHandler::add_radiation_contribution for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { long const np = pti.numParticles(); auto& attribs = pti.GetAttribs(); - auto& p_w = attribs[PIdx::w]; + auto& p_w = attribs[PIdx::w]; auto& p_ux = attribs[PIdx::ux]; auto& p_uy = attribs[PIdx::uy]; auto& p_uz = attribs[PIdx::uz]; - auto& attribut = attribs[PIdx::nattribs]; - + auto& attribut = attribs[PIdx::nattribs+1]; + amrex::Print() << attribut[-1] << "\n"; //auto& ux = src.m_rdata[PIdx::ux][i_src]; //auto& uy = src.m_rdata[PIdx::uy][i_src]; //auto& uz = src.m_rdata[PIdx::uz][i_src]; @@ -110,16 +109,16 @@ void RadiationHandler::add_detector // Box of angle and frequency const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); - amrex::Print() << detect_box << "\n"; - ncomp = 0; + ncomp = 2; //Real and complex part amrex::FArrayBox fab_detect(detect_box, ncomp); - amrex::Print() << fab_detect << "\n"; //Calculation of angle resolution + m_d_theta.resize(2); for(int i=0; i<2; i++){ m_d_theta[i] = 2*m_theta_range[i]/static_cast(m_det_pts[i]); } //Set the resolution of the detector + m_d_d.resize(2); for(int idim = 0; idim<3; ++idim){ if(m_det_direction[idim]==1){ m_d_d[(idim+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); From abdde489ddaf184c61bf6a5aebb41b303c40f225 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Wed, 6 Dec 2023 15:26:37 +0900 Subject: [PATCH 07/60] . --- Source/Evolve/WarpXEvolve.cpp | 2 +- Source/Particles/MultiParticleContainer.cpp | 4 +- .../NamedComponentParticleContainer.H | 2 +- .../Particles/Radiation/RadiationHandler.cpp | 75 +++++++++++++++---- 4 files changed, 63 insertions(+), 20 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index e7e910cef59..068db696914 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -432,7 +432,7 @@ WarpX::OneStep_nosub (Real cur_time) //Radiation contribution at each timestep //Only level 0 is supported - mypc->doRadiation(dt[0]); + mypc->doRadiation(dt[0],cur_time); diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 3bcf4a9781c..2b3cf2d843e 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -960,8 +960,8 @@ void MultiParticleContainer::keepoldmomentum(){ //Load the momentums amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); - amrex::ParticleReal* uy = soa.GetRealData(PIdx::ux).data(); - amrex::ParticleReal* uz = soa.GetRealData(PIdx::ux).data(); + amrex::ParticleReal* uy = soa.GetRealData(PIdx::uy).data(); + amrex::ParticleReal* uz = soa.GetRealData(PIdx::uz).data(); //Finding the good attribute index int index_name_x=pc->GetRealCompIndex("prev_u_x"); diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index 7986cb0cc31..8e984407aa3 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -183,7 +183,7 @@ public: */ int GetRealCompIndex (const std::string& name) { - return particle_runtime_comps.at(name); + return particle_comps.at(name); } /** Initialize value at a component * diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index ad511ceb674..e3aa6a94558 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -11,6 +11,9 @@ #include +#include + + #include using namespace warpx::utils; @@ -51,10 +54,8 @@ RadiationHandler::RadiationHandler() } void RadiationHandler::add_radiation_contribution - (const amrex::Real dt, std::unique_ptr& pc){ - const auto level0=0; - const auto part=pc->GetParticles(level0); - + (const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time){ + const auto level0=0; #ifdef AMREX_USE_OMP #pragma omp parallel for #endif @@ -66,24 +67,66 @@ void RadiationHandler::add_radiation_contribution auto& p_ux = attribs[PIdx::ux]; auto& p_uy = attribs[PIdx::uy]; auto& p_uz = attribs[PIdx::uz]; - auto& attribut = attribs[PIdx::nattribs+1]; - amrex::Print() << attribut[-1] << "\n"; - //auto& ux = src.m_rdata[PIdx::ux][i_src]; - //auto& uy = src.m_rdata[PIdx::uy][i_src]; - //auto& uz = src.m_rdata[PIdx::uz][i_src]; - //auto GetPosition = GetParticlePosition(pti); + auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + auto& part=pc->GetParticles(level0)[index]; + auto& soa = part.GetStructOfArrays(); + + const auto prev_u_x_idx =pc->GetRealCompIndex("prev_u_x"); + const auto prev_u_y_idx =pc->GetRealCompIndex("prev_u_y"); + const auto prev_u_z_idx =pc->GetRealCompIndex("prev_u_z"); + + auto* p_ux_old = soa.GetRealData(prev_u_x_idx).data(); + auto* p_uy_old = soa.GetRealData(prev_u_y_idx).data(); + auto* p_uz_old = soa.GetRealData(prev_u_z_idx).data(); + + auto GetPosition = GetParticlePosition(pti); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip) - { //amrex::ParticleReal xp, yp, zp; - //GetPosition.AsStored(ip, xp, yp, zp); - const amrex::ParticleReal p_ux_unit=p_ux[ip]; - const amrex::ParticleReal p_uy_unit=p_uy[ip]; - const amrex::ParticleReal p_uz_unit=p_uz[ip]; + { amrex::ParticleReal xp, yp, zp; + GetPosition.AsStored(ip, xp, yp, zp); + + std::cout << "######### " << ip << " " + << xp << " " << yp << " " << zp << " " + << p_ux[ip] << " " << p_uy[ip] << " " << p_uz[ip] << " " + << p_ux_old[ip] << " " << p_uy_old[ip] << " " << p_uz_old[ip] << std::endl; + + amrex::ParticleReal p_ux_unit=p_ux_old[ip]; + amrex::ParticleReal p_uy_unit=p_uy_old[ip]; + amrex::ParticleReal p_uz_unit=p_uz_old[ip]; + + amrex::ParticleReal p_ux_old_unit=p_ux[ip]; + amrex::ParticleReal p_uy_old_unit=p_uy[ip]; + amrex::ParticleReal p_uz_old_unit=p_uz[ip]; amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); + //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal + amrex::ParticleReal un_betan = 1-(p_ux_unit*static_cast(m_det_direction[0])+p_uy_unit*static_cast(m_det_direction[1])+p_uz_unit*static_cast(m_det_direction[2])); + + //Calculation of gamma + amrex::ParticleReal gamma = 1/(1-std::pow(p_u,2)/std::pow(ablastr::constant::SI::c,2)); + + //Calculation of betapoint + amrex::ParticleReal betapointx=(p_ux_unit-p_ux_old_unit)/dt/ablastr::constant::SI::c; + amrex::ParticleReal betapointy=(p_uy_unit-p_uy_old_unit)/dt/ablastr::constant::SI::c; + amrex::ParticleReal betapointz=(p_uz_unit-p_uz_old_unit)/dt/ablastr::constant::SI::c; + + //Calculation of nxbeta + amrex::ParticleReal ncrossBetapointx=(static_cast(m_det_direction[1])-p_uy_unit/gamma)*betapointz - (static_cast(m_det_direction[2])-p_uz_unit/gamma)*betapointy; + amrex::ParticleReal ncrossBetapointy=(static_cast(m_det_direction[2])-p_uz_unit/gamma)*betapointx - (static_cast(m_det_direction[0])-p_ux_unit/gamma)*betapointz; + amrex::ParticleReal ncrossBetapointz=(static_cast(m_det_direction[0])-p_ux_unit/gamma)*betapointy - (static_cast(m_det_direction[1])-p_uy_unit/gamma)*betapointx; + + //Calculation of nxnxbeta + amrex::ParticleReal ncrossncrossBetapointx=static_cast(m_det_direction[2])*ncrossBetapointz-static_cast(m_det_direction[2])*ncrossBetapointy; + amrex::ParticleReal ncrossncrossBetapointy=static_cast(m_det_direction[0])*ncrossBetapointx-static_cast(m_det_direction[0])*ncrossBetapointz; + amrex::ParticleReal ncrossncrossBetapointz=static_cast(m_det_direction[1])*ncrossBetapointy-static_cast(m_det_direction[1])*ncrossBetapointx; + + //Calculation of ei(omegat-n.r) + for(int i_om=0, i_om Date: Wed, 6 Dec 2023 15:50:27 +0900 Subject: [PATCH 08/60] . --- Source/Particles/Radiation/RadiationHandler.cpp | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index e3aa6a94558..6694c084521 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -122,10 +122,17 @@ void RadiationHandler::add_radiation_contribution amrex::ParticleReal ncrossncrossBetapointz=static_cast(m_det_direction[1])*ncrossBetapointy-static_cast(m_det_direction[1])*ncrossBetapointx; //Calculation of ei(omegat-n.r) + + //omega effective calculé + amrex::Real omega_calc = 0 for(int i_om=0, i_om Date: Thu, 7 Dec 2023 09:07:40 +0900 Subject: [PATCH 09/60] . --- Source/Particles/MultiParticleContainer.cpp | 1 - Source/Particles/Radiation/RadiationHandler.H | 2 + .../Particles/Radiation/RadiationHandler.cpp | 73 ++++++++++++------- 3 files changed, 49 insertions(+), 27 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 2b3cf2d843e..8a8ff313ecf 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1031,7 +1031,6 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt) if (pc->has_radiation()){ m_p_radiation_handler->add_radiation_contribution(dt,pc); - } } } diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index ff05811625b..70c3d936d09 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -33,6 +33,8 @@ public: //void doRadiation (amrex::Real dt, MultiParticleContainer* mypc); void add_radiation_contribution(const amrex::Real dt, std::unique_ptr& pc); void add_detector(); + void RadiationHandler::gather_and_write_radiation(const std::string& filename); + private: diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 6694c084521..7330c12892c 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -91,13 +91,15 @@ void RadiationHandler::add_radiation_contribution << p_ux[ip] << " " << p_uy[ip] << " " << p_uz[ip] << " " << p_ux_old[ip] << " " << p_uy_old[ip] << " " << p_uz_old[ip] << std::endl; - amrex::ParticleReal p_ux_unit=p_ux_old[ip]; - amrex::ParticleReal p_uy_unit=p_uy_old[ip]; - amrex::ParticleReal p_uz_unit=p_uz_old[ip]; + amrex::ParticleReal p_ux_old_unit=p_ux_old[ip]; + amrex::ParticleReal p_uy_old_unit=p_uy_old[ip]; + amrex::ParticleReal p_uz_old_unit=p_uz_old[ip]; - amrex::ParticleReal p_ux_old_unit=p_ux[ip]; - amrex::ParticleReal p_uy_old_unit=p_uy[ip]; - amrex::ParticleReal p_uz_old_unit=p_uz[ip]; + amrex::ParticleReal p_ux_unit=p_ux[ip]; + amrex::ParticleReal p_uy_unit=p_uy[ip]; + amrex::ParticleReal p_uz_unit=p_uz[ip]; + + amrex::ParticleReal const q = pc.getCharge(); amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal @@ -129,31 +131,50 @@ void RadiationHandler::add_radiation_contribution for(int i_x=0, i_x(m_det_direction[(idim+1)%3])*(xp-m_det_distance*static_cast(m_det_direction[0])-pos_det_x[0])-static_cast(m_det_direction[1])*(yp-m_det_distance*static_cast(m_det_direction[1]))-static_cast(m_det_direction[2])*(zp-m_det_distance*static_cast(m_det_direction[2])))/ablastr::constant::SI::c); + m_d_d[(idim+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); + m_d_d[(idim+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); + } + amrex::Complex eiomega=(amrex::cos(dephas), amrex::sin(dephas)); + amrex::Complex Term_x=ncrossncrossBetapointx/un_betan**2*eiomega; + amrex::Complex Term_y=ncrossncrossBetapointy/un_betan**2*eiomega; + amrex::Complex Term_z=ncrossncrossBetapointz/un_betan**2*eiomega; + + //Calcul du terme direct de contribution + amrex::Real Term_int = q*dt/(16*ablastr::constant::SI::pi*ablastr::constant::SI::mu0*ablastr::constant::SI::cx)*std(std::pow(qpTerm_x,2)+std::pow(Term_y,2)+std::pow(Term_z,2)); + fab_detect.saxpy(Real 1, fab_detect); + gather_and_write_radiation(const std::string& filename) + } + } } + } + } + } + + }); + +void RadiationHandler::gather_and_write_radiation(const std::string& filename) +{ + auto radiation_data_cpu = amrex::Vector(m_I*m_J*m_W*2); + amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, + m_radiation_data.begin(), m_radiation_data.end(), radiation_data_cpu.begin()); + amrex::Gpu::streamSynchronize(); + amrex::ParallelDescriptor::ReduceRealSum(radiation_data_cpu.data(), radiation_data_cpu.size()); - }); - -// const int total_partdiag_size = amrex::Scan::ExclusiveSum(np,Flag,IndexLocation); -// auto& ptile_dst = pc_dst.DefineAndReturnParticleTile(lev, pti.index(), pti.LocalTileIndex() ); -// auto old_size = ptile_dst.numParticles(); -// ptile_dst.resize(old_size + total_partdiag_size); -// amrex::filterParticles(ptile_dst, ptile_src, GetParticleFilter, 0, old_size, np); -// auto dst_data = ptile_dst.getParticleTileData(); -// amrex::ParallelFor(np, -// [=] AMREX_GPU_DEVICE(int i) -// { -// if (Flag[i] == 1) GetParticleLorentzTransform(dst_data, src_data, i, -// old_size + IndexLocation[i]); -// }); -// amrex::Gpu::synchronize(); + if (amrex::ParallelDescriptor::IOProcessor() ){ + auto of = std::ofstream(filename, std::ios::binary); + + for (const auto& d : radiation_data_cpu){ + of << static_cast(d); } + + of.close(); } } - void RadiationHandler::add_detector (){ @@ -161,7 +182,7 @@ void RadiationHandler::add_detector const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); ncomp = 2; //Real and complex part amrex::FArrayBox fab_detect(detect_box, ncomp); - + fab_detect.setval(0,0) //Calculation of angle resolution m_d_theta.resize(2); for(int i=0; i<2; i++){ From c3bf08a92395bc1c71aba629f47f2875cc0edfbf Mon Sep 17 00:00:00 2001 From: tmsclark Date: Thu, 7 Dec 2023 13:52:02 +0900 Subject: [PATCH 10/60] . --- Source/Evolve/WarpXEvolve.cpp | 2 + Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 4 +- Source/Particles/Radiation/RadiationHandler.H | 18 +++- .../Particles/Radiation/RadiationHandler.cpp | 99 +++++++++++-------- Source/Particles/Sorting/SortingUtils.H | 9 ++ Source/Particles/Sorting/SortingUtils.cpp | 12 +++ 7 files changed, 97 insertions(+), 49 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 068db696914..0927a712804 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 @@ -244,6 +245,7 @@ WarpX::Evolve (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1], *Bfield_aux[lev][2]); } + //gather_and_write_radiation( filename); is_synchronized = true; } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index a931a00f636..7243c4e606d 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -330,7 +330,7 @@ public: const amrex::MultiFab& Bz); #endif - void doRadiation(amrex::Real dt); + void doRadiation(amrex::Real dt, amrex::Real cur_time); void keepoldmomentum(); void keepoldmomentum_p(std::unique_ptr& pc); diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 8a8ff313ecf..26361e0245a 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1024,13 +1024,13 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vectorhas_radiation()){ - m_p_radiation_handler->add_radiation_contribution(dt,pc); + m_p_radiation_handler->add_radiation_contribution(dt,pc,cur_time); } } } diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 70c3d936d09..68be8055dae 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -17,7 +17,7 @@ #include #include #include - +#include "Particles/Sorting/SortingUtils.H" #include #include @@ -31,9 +31,9 @@ public: /* 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); - void add_detector(); - void RadiationHandler::gather_and_write_radiation(const std::string& filename); + void add_radiation_contribution(const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time); + amrex::IntVect add_detector(); + void gather_and_write_radiation(const std::string& filename); private: @@ -65,6 +65,16 @@ private: //resolution of the grid of the detectir amrex::Vector m_d_d; + // + amrex::Gpu::DeviceVector m_radiation_data; + + // + amrex::Real dephas; + + // + amrex::Vector> det_pos; + amrex::Vector omega_calc; + }; #endif // WARPX_PARTICLES_RADIATION_H_ diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 7330c12892c..444c78f5adf 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -42,7 +42,22 @@ RadiationHandler::RadiationHandler() pp_radiations.getarr("detector_number_points", m_det_pts,0,2); pp_radiations.query("detector_direction", m_det_direction); pp_radiations.query("detector_distance", m_det_distance); - add_detector(); + amrex::IntVect d_d_fab = add_detector(); + /* Initialize the Fab with the field in function of angle and frequency */ + //int numcomps = 4; + //BaseFab fab(bx,numcomps); + // Box of angle and frequency + //const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); + //ncomp = 2; //Real and complex part + //amrex::FArrayBox fab_detect(detect_box, ncomp); + //fab_detect.setval(0,0) + m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); + fillWithConsecutiveReal(det_pos[0], 0,s); + fillWithConsecutiveReal(det_pos[1]); + fillWithConsecutiveReal(omega_calc); + + + } // @@ -81,15 +96,14 @@ void RadiationHandler::add_radiation_contribution auto* p_uz_old = soa.GetRealData(prev_u_z_idx).data(); auto GetPosition = GetParticlePosition(pti); + amrex::ParticleReal const q = pc->getCharge(); + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip) { amrex::ParticleReal xp, yp, zp; - GetPosition.AsStored(ip, xp, yp, zp); - - std::cout << "######### " << ip << " " - << xp << " " << yp << " " << zp << " " - << p_ux[ip] << " " << p_uy[ip] << " " << p_uz[ip] << " " - << p_ux_old[ip] << " " << p_uy_old[ip] << " " << p_uz_old[ip] << std::endl; + + GetPosition.AsStored(ip,xp, yp, zp); + amrex::GpuArray Part_pos{xp,yp,zp}; amrex::ParticleReal p_ux_old_unit=p_ux_old[ip]; amrex::ParticleReal p_uy_old_unit=p_uy_old[ip]; @@ -99,8 +113,6 @@ void RadiationHandler::add_radiation_contribution amrex::ParticleReal p_uy_unit=p_uy[ip]; amrex::ParticleReal p_uz_unit=p_uz[ip]; - amrex::ParticleReal const q = pc.getCharge(); - amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal amrex::ParticleReal un_betan = 1-(p_ux_unit*static_cast(m_det_direction[0])+p_uy_unit*static_cast(m_det_direction[1])+p_uz_unit*static_cast(m_det_direction[2])); @@ -126,39 +138,37 @@ void RadiationHandler::add_radiation_contribution //Calculation of ei(omegat-n.r) //omega effective calculé - amrex::Real omega_calc = 0 - for(int i_om=0, i_om(m_det_direction[(idim+1)%3])*(xp-m_det_distance*static_cast(m_det_direction[0])-pos_det_x[0])-static_cast(m_det_direction[1])*(yp-m_det_distance*static_cast(m_det_direction[1]))-static_cast(m_det_direction[2])*(zp-m_det_distance*static_cast(m_det_direction[2])))/ablastr::constant::SI::c); - m_d_d[(idim+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); - m_d_d[(idim+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); + amrex::Real omega_calc = 0; + for(int i_om=0; i_om(m_I*m_J*m_W*2); + auto radiation_data_cpu = amrex::Vector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, m_radiation_data.begin(), m_radiation_data.end(), radiation_data_cpu.begin()); amrex::Gpu::streamSynchronize(); @@ -175,14 +185,9 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) of.close(); } } -void RadiationHandler::add_detector +amrex::IntVect RadiationHandler::add_detector (){ - - // Box of angle and frequency - const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); - ncomp = 2; //Real and complex part - amrex::FArrayBox fab_detect(detect_box, ncomp); - fab_detect.setval(0,0) + Geometry const& geom = WarpX::GetInstance().Geom(lev); //Calculation of angle resolution m_d_theta.resize(2); for(int i=0; i<2; i++){ @@ -190,6 +195,7 @@ void RadiationHandler::add_detector } //Set the resolution of the detector m_d_d.resize(2); + geom[0].ProbLo(); for(int idim = 0; idim<3; ++idim){ if(m_det_direction[idim]==1){ m_d_d[(idim+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); @@ -197,5 +203,14 @@ void RadiationHandler::add_detector } m_d_omega=m_omega_range[1]-m_omega_range[0]/static_cast(m_omega_points); + + + //Calculate the sides of the detector + size_dim = + //fillWithConsecutiveReal(pos_det_x) + //pos_det_y + //omega_calc + } + return(amrex::IntVect(m_d_d[0],m_d_d[1],m_d_omega)); } \ No newline at end of file diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index ac2c63e88f8..076c6a21e0c 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -21,6 +21,15 @@ */ void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ); +/** \brief Fill the elements of the input vector with consecutive integer, + * starting from 0 + * + * \param[inout] v Vector of reals, to be filled by this routine + */ + +void fillWithConsecutiveReal( amrex::Gpu::DeviceVector& 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/Sorting/SortingUtils.cpp b/Source/Particles/Sorting/SortingUtils.cpp index 699119e8e18..746967604e1 100644 --- a/Source/Particles/Sorting/SortingUtils.cpp +++ b/Source/Particles/Sorting/SortingUtils.cpp @@ -20,3 +20,15 @@ void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ) std::iota( v.begin(), v.end(), 0L ); #endif } + +void fillWithConsecutiveReal( amrex::Gpu::DeviceVector& v,amrex::Real begin,amrex::Real increment, int N) +{ +#ifdef AMREX_USE_GPU + // On GPU: Use amrex + auto data = begin; + AMREX_FOR_1D( N, i, {data+i*increment;}); +#else + // On CPU: Use std library + std::iota( begin, N, 0L ); +#endif +} From 7db336a5df29022dda0261e5f14db996f31301b2 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Fri, 15 Dec 2023 13:49:45 +0100 Subject: [PATCH 11/60] . --- Source/Evolve/WarpXEvolve.cpp | 2 +- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 8 + Source/Particles/Radiation/RadiationHandler.H | 13 +- .../Particles/Radiation/RadiationHandler.cpp | 171 +++++++++++------- Source/Particles/Sorting/SortingUtils.cpp | 12 +- 6 files changed, 124 insertions(+), 84 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 0927a712804..b6cdaec8139 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -245,7 +245,7 @@ WarpX::Evolve (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1], *Bfield_aux[lev][2]); } - //gather_and_write_radiation( filename); + mypc->Dump_radiations(); is_synchronized = true; } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 7243c4e606d..6d7cd0d6a8d 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -335,7 +335,7 @@ public: void keepoldmomentum(); void keepoldmomentum_p(std::unique_ptr& pc); - + void Dump_radiations(); int getSpeciesID (std::string product_str) const; protected: diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 26361e0245a..6b152c2e1d7 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1035,6 +1035,14 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt, amrex::Real cur_ } } } + + +void MultiParticleContainer::Dump_radiations(){ + if (m_at_least_one_has_radiation){ + m_p_radiation_handler->Integral_overtime(); + m_p_radiation_handler->gather_and_write_radiation("Radiation"); + } +} #ifdef WARPX_QED void MultiParticleContainer::InitQED () { diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 68be8055dae..b4117516f65 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -17,6 +17,9 @@ #include #include #include +#include + +#include "WarpX.H" #include "Particles/Sorting/SortingUtils.H" #include #include @@ -32,8 +35,9 @@ public: /* 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); - amrex::IntVect add_detector(); + void add_detector(); void gather_and_write_radiation(const std::string& filename); + void Integral_overtime(); private: @@ -53,25 +57,28 @@ private: double m_d_omega; //The vector with the resolution of the detector amrex::Vector m_d_det; - int m_det_distance; + amrex::Real m_det_distance; //Define the Fab with the datas int ncomp; amrex::Vector m_det_pts; - amrex::IntVect m_det_direction; + amrex::Vector m_det_direction; amrex::Vector m_pos_det; + amrex::Vector n; //resolution of the grid of the detectir amrex::Vector m_d_d; // amrex::Gpu::DeviceVector m_radiation_data; + amrex::Gpu::DeviceVector m_radiation_calculation; // amrex::Real dephas; // + amrex::Vector> det_bornes; amrex::Vector> det_pos; amrex::Vector omega_calc; diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 444c78f5adf..8f96e7d4716 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -40,39 +40,44 @@ RadiationHandler::RadiationHandler() //Type of detector pp_radiations.getarr("detector_number_points", m_det_pts,0,2); - pp_radiations.query("detector_direction", m_det_direction); + pp_radiations.getarr("detector_direction", m_det_direction,0,3); pp_radiations.query("detector_distance", m_det_distance); - amrex::IntVect d_d_fab = add_detector(); + add_detector(); /* Initialize the Fab with the field in function of angle and frequency */ //int numcomps = 4; //BaseFab fab(bx,numcomps); // Box of angle and frequency //const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); - //ncomp = 2; //Real and complex part + ncomp = 6; //Real and complex part //amrex::FArrayBox fab_detect(detect_box, ncomp); //fab_detect.setval(0,0) m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); - fillWithConsecutiveReal(det_pos[0], 0,s); - fillWithConsecutiveReal(det_pos[1]); - fillWithConsecutiveReal(omega_calc); - + m_radiation_calculation = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points); + det_pos.resize(2); + det_pos[0].resize(m_det_pts[0]); + det_pos[1].resize(m_det_pts[1]); + omega_calc.resize(m_omega_points); - } - - // - + //Initialization : + for(int i=0; i fab(bx,numcomps); + } } void RadiationHandler::add_radiation_contribution (const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time){ const auto level0=0; #ifdef AMREX_USE_OMP -#pragma omp parallel for +#pragma omp parallel #endif { for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { @@ -101,7 +106,7 @@ void RadiationHandler::add_radiation_contribution amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip) { amrex::ParticleReal xp, yp, zp; - + GetPosition.AsStored(ip,xp, yp, zp); amrex::GpuArray Part_pos{xp,yp,zp}; @@ -112,65 +117,74 @@ void RadiationHandler::add_radiation_contribution amrex::ParticleReal p_ux_unit=p_ux[ip]; amrex::ParticleReal p_uy_unit=p_uy[ip]; amrex::ParticleReal p_uz_unit=p_uz[ip]; - + amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); - //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal - amrex::ParticleReal un_betan = 1-(p_ux_unit*static_cast(m_det_direction[0])+p_uy_unit*static_cast(m_det_direction[1])+p_uz_unit*static_cast(m_det_direction[2])); //Calculation of gamma amrex::ParticleReal gamma = 1/(1-std::pow(p_u,2)/std::pow(ablastr::constant::SI::c,2)); - //Calculation of betapoint - amrex::ParticleReal betapointx=(p_ux_unit-p_ux_old_unit)/dt/ablastr::constant::SI::c; - amrex::ParticleReal betapointy=(p_uy_unit-p_uy_old_unit)/dt/ablastr::constant::SI::c; - amrex::ParticleReal betapointz=(p_uz_unit-p_uz_old_unit)/dt/ablastr::constant::SI::c; - - //Calculation of nxbeta - amrex::ParticleReal ncrossBetapointx=(static_cast(m_det_direction[1])-p_uy_unit/gamma)*betapointz - (static_cast(m_det_direction[2])-p_uz_unit/gamma)*betapointy; - amrex::ParticleReal ncrossBetapointy=(static_cast(m_det_direction[2])-p_uz_unit/gamma)*betapointx - (static_cast(m_det_direction[0])-p_ux_unit/gamma)*betapointz; - amrex::ParticleReal ncrossBetapointz=(static_cast(m_det_direction[0])-p_ux_unit/gamma)*betapointy - (static_cast(m_det_direction[1])-p_uy_unit/gamma)*betapointx; - - //Calculation of nxnxbeta - amrex::ParticleReal ncrossncrossBetapointx=static_cast(m_det_direction[2])*ncrossBetapointz-static_cast(m_det_direction[2])*ncrossBetapointy; - amrex::ParticleReal ncrossncrossBetapointy=static_cast(m_det_direction[0])*ncrossBetapointx-static_cast(m_det_direction[0])*ncrossBetapointz; - amrex::ParticleReal ncrossncrossBetapointz=static_cast(m_det_direction[1])*ncrossBetapointy-static_cast(m_det_direction[1])*ncrossBetapointx; - //Calculation of ei(omegat-n.r) - + n.resize(3); //omega effective calculé - amrex::Real omega_calc = 0; for(int i_om=0; i_om(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); + 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_data.begin(), m_radiation_data.end(), radiation_data_cpu.begin()); + 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()); @@ -185,32 +199,53 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) of.close(); } } -amrex::IntVect RadiationHandler::add_detector +void RadiationHandler::add_detector (){ - Geometry const& geom = WarpX::GetInstance().Geom(lev); + const auto level0=0; + amrex::Geometry const& geom = WarpX::GetInstance().Geom(level0); //Calculation of angle resolution - m_d_theta.resize(2); + m_d_theta.resize(2); for(int i=0; i<2; i++){ m_d_theta[i] = 2*m_theta_range[i]/static_cast(m_det_pts[i]); } //Set the resolution of the detector - m_d_d.resize(2); - geom[0].ProbLo(); - for(int idim = 0; idim<3; ++idim){ - if(m_det_direction[idim]==1){ - m_d_d[(idim+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); - m_d_d[(idim+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); + m_d_d.resize(3); + amrex::Vector center; + center.resize(3); + for(int idim=0; idim(m_omega_points); - //Calculate the sides of the detector - size_dim = + det_bornes.resize(2); + det_bornes[0].resize(2); + det_bornes[1].resize(2); + + for(int idim = 0; idim<2; ++idim){ + det_bornes[idim][0]=center[idim]-std::tan(m_theta_range[idim]/2)*m_det_distance; + det_bornes[idim][1]=center[idim]+std::tan(m_theta_range[idim]/2)*m_det_distance; + } //fillWithConsecutiveReal(pos_det_x) //pos_det_y //omega_calc } - return(amrex::IntVect(m_d_d[0],m_d_d[1],m_d_omega)); +} + +void RadiationHandler::Integral_overtime() +{ + for(int i_om=0; i_om& v ) #endif } -void fillWithConsecutiveReal( amrex::Gpu::DeviceVector& v,amrex::Real begin,amrex::Real increment, int N) -{ -#ifdef AMREX_USE_GPU - // On GPU: Use amrex - auto data = begin; - AMREX_FOR_1D( N, i, {data+i*increment;}); -#else - // On CPU: Use std library - std::iota( begin, N, 0L ); -#endif -} + From c84b257f5896e3c44e255f8292b3477c1da2ec06 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Fri, 15 Dec 2023 14:40:40 +0100 Subject: [PATCH 12/60] . --- Source/Particles/Radiation/RadiationHandler.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 8f96e7d4716..abc64bafe61 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -135,8 +135,8 @@ void RadiationHandler::add_radiation_contribution } else{ - n[idimo]=(Part_pos[idimo]-det_pos[(idimo+1)%3][i_y])/m_det_distance; - n[idimo]=(Part_pos[idimo]-det_pos[(idimo+2)%3][i_x])/m_det_distance; + n[idimo]=(Part_pos[idimo]-det_pos[(idimo+1)%2][i_y])/m_det_distance; + n[idimo]=(Part_pos[idimo]-det_pos[(idimo+2)%2][i_x])/m_det_distance; } } //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal From b8a2aa8b2286bae7cc8aa711d3e963d9e6ce374d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 15 Dec 2023 13:47:31 +0000 Subject: [PATCH 13/60] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- Source/Evolve/WarpXEvolve.cpp | 6 ++-- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 8 ++--- .../NamedComponentParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/Radiation/RadiationHandler.H | 11 +++---- .../Particles/Radiation/RadiationHandler.cpp | 32 +++++++++---------- Source/Particles/Sorting/SortingUtils.cpp | 2 -- Source/Utils/Math/Complex.H | 14 ++++---- Source/WarpX.H | 2 +- Source/WarpX.cpp | 2 +- 11 files changed, 40 insertions(+), 43 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index d89de935108..2d5101d7cdc 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -434,13 +434,13 @@ WarpX::OneStep_nosub (Real cur_time) PushParticlesandDeposit(cur_time); //Radiation contribution at each timestep - //Only level 0 is supported + //Only level 0 is supported mypc->doRadiation(dt[0],cur_time); - - + + ExecutePythonCallback("afterdeposition"); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 6b13415c5b1..016b4e84380 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -488,7 +488,7 @@ protected: amrex::Real m_qed_schwinger_zmax = std::numeric_limits::max(); /** - * For Radiation module + * For Radiation module */ bool m_at_least_one_particle_radiate = 0; int m_nspecies_radiate; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index efd3705f271..aaf3bcaf315 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -377,7 +377,7 @@ MultiParticleContainer::ReadParameters () m_qed_schwinger_threshold_poisson_gaussian); utils::parser::queryWithParser( pp_qed_schwinger, "xmin", m_qed_schwinger_xmin); - + #if defined(WARPX_DIM_3D) utils::parser::queryWithParser( pp_qed_schwinger, "ymin", m_qed_schwinger_ymin); @@ -962,7 +962,7 @@ void MultiParticleContainer::keepoldmomentum(){ }); } - } + } void MultiParticleContainer::doCollisions ( Real cur_time, amrex::Real dt ) @@ -1005,11 +1005,11 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vectorhas_radiation()){ - + m_p_radiation_handler->add_radiation_contribution(dt,pc,cur_time); } } diff --git a/Source/Particles/NamedComponentParticleContainer.H b/Source/Particles/NamedComponentParticleContainer.H index 924a493d899..394b808ad8b 100644 --- a/Source/Particles/NamedComponentParticleContainer.H +++ b/Source/Particles/NamedComponentParticleContainer.H @@ -176,7 +176,7 @@ public: [[nodiscard]] std::map getParticleRuntimeComps () const noexcept { return particle_runtime_comps;} /** 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 diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 356e7722498..0453a626071 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -359,7 +359,7 @@ public: #endif bool has_radiation() const override {return m_do_radiation;} - + protected: std::string species_name; std::vector> plasma_injectors; diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index b4117516f65..4d195005879 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -33,7 +33,7 @@ public: RadiationHandler (); /* Perform the calculation of the radiation */ - //void doRadiation (amrex::Real dt, MultiParticleContainer* mypc); + //void doRadiation (amrex::Real dt, MultiParticleContainer* mypc); void add_radiation_contribution(const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time); void add_detector(); void gather_and_write_radiation(const std::string& filename); @@ -44,13 +44,13 @@ private: // Frequency range amrex::Vector m_omega_range; - - // Dimensions of the detector + + // Dimensions of the detector amrex::Vector m_theta_range; int m_omega_points; - //put a detector + //put a detector int m_get_a_detector; // Resolution of the detector amrex::Vector m_d_theta; @@ -67,7 +67,7 @@ private: amrex::Vector m_pos_det; amrex::Vector n; - //resolution of the grid of the detectir + //resolution of the grid of the detectir amrex::Vector m_d_d; // @@ -84,4 +84,3 @@ private: }; #endif // WARPX_PARTICLES_RADIATION_H_ - diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index abc64bafe61..02d58e0e40e 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -18,17 +18,17 @@ using namespace warpx::utils; RadiationHandler::RadiationHandler() -{ - // Read in radiation input +{ + // Read in radiation input const amrex::ParmParse pp_radiations("radiations"); - // Verify if there is a detector + // Verify if there is a detector pp_radiations.query("put_a_detector", m_get_a_detector); #if defined WARPX_DIM_RZ WARPX_ABORT_WITH_MESSAGE("Radiation is not supported yet with RZ."); -#endif +#endif //Compute the radiation if(m_get_a_detector){ @@ -75,9 +75,9 @@ RadiationHandler::RadiationHandler() void RadiationHandler::add_radiation_contribution (const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time){ - const auto level0=0; + const auto level0=0; #ifdef AMREX_USE_OMP -#pragma omp parallel +#pragma omp parallel #endif { for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { @@ -117,10 +117,10 @@ void RadiationHandler::add_radiation_contribution amrex::ParticleReal p_ux_unit=p_ux[ip]; amrex::ParticleReal p_uy_unit=p_uy[ip]; amrex::ParticleReal p_uz_unit=p_uz[ip]; - + amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); - //Calculation of gamma + //Calculation of gamma amrex::ParticleReal gamma = 1/(1-std::pow(p_u,2)/std::pow(ablastr::constant::SI::c,2)); //Calculation of ei(omegat-n.r) @@ -136,7 +136,7 @@ void RadiationHandler::add_radiation_contribution } else{ n[idimo]=(Part_pos[idimo]-det_pos[(idimo+1)%2][i_y])/m_det_distance; - n[idimo]=(Part_pos[idimo]-det_pos[(idimo+2)%2][i_x])/m_det_distance; + n[idimo]=(Part_pos[idimo]-det_pos[(idimo+2)%2][i_x])/m_det_distance; } } //Calculation of 1_beta.n, n corresponds to m_det_direction, the direction of the normal @@ -156,7 +156,7 @@ void RadiationHandler::add_radiation_contribution amrex::ParticleReal ncrossncrossBetapointx=n[2]*ncrossBetapointz-n[2]*ncrossBetapointy; amrex::ParticleReal ncrossncrossBetapointy=n[0]*ncrossBetapointx-n[0]*ncrossBetapointz; amrex::ParticleReal ncrossncrossBetapointz=n[1]*ncrossBetapointy-n[1]*ncrossBetapointx; - + //function cospi and not cos Complex eiomega=(amrex::Math::cospi(dephas/ablastr::constant::math::pi), amrex::Math::cospi((ablastr::constant::math::pi/2-dephas)/ablastr::constant::math::pi)); Complex Term_x= q*dt/(16*pow(ablastr::constant::math::pi,3)*ablastr::constant::SI::ep0*ablastr::constant::SI::c)*ncrossncrossBetapointx/std::pow(un_betan,2)*eiomega; @@ -174,11 +174,11 @@ void RadiationHandler::add_radiation_contribution } } } - + }); } } - } + } void RadiationHandler::gather_and_write_radiation(const std::string& filename) { @@ -203,7 +203,7 @@ void RadiationHandler::add_detector (){ const auto level0=0; amrex::Geometry const& geom = WarpX::GetInstance().Geom(level0); - //Calculation of angle resolution + //Calculation of angle resolution m_d_theta.resize(2); for(int i=0; i<2; i++){ m_d_theta[i] = 2*m_theta_range[i]/static_cast(m_det_pts[i]); @@ -223,7 +223,7 @@ void RadiationHandler::add_detector m_d_omega=m_omega_range[1]-m_omega_range[0]/static_cast(m_omega_points); - //Calculate the sides of the detector + //Calculate the sides of the detector det_bornes.resize(2); det_bornes[0].resize(2); det_bornes[1].resize(2); @@ -235,7 +235,7 @@ void RadiationHandler::add_detector //fillWithConsecutiveReal(pos_det_x) //pos_det_y //omega_calc - + } } @@ -248,4 +248,4 @@ void RadiationHandler::Integral_overtime() } } } -} \ No newline at end of file +} diff --git a/Source/Particles/Sorting/SortingUtils.cpp b/Source/Particles/Sorting/SortingUtils.cpp index 9880a2d7fbc..699119e8e18 100644 --- a/Source/Particles/Sorting/SortingUtils.cpp +++ b/Source/Particles/Sorting/SortingUtils.cpp @@ -20,5 +20,3 @@ void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ) std::iota( v.begin(), v.end(), 0L ); #endif } - - diff --git a/Source/Utils/Math/Complex.H b/Source/Utils/Math/Complex.H index 81a29407a31..651dc3449dc 100644 --- a/Source/Utils/Math/Complex.H +++ b/Source/Utils/Math/Complex.H @@ -14,15 +14,15 @@ #if defined(AMREX_USE_CUDA) # include #elif defined(AMREX_USE_HIP) -# include +# include #else # if WARPX_USE_PSATD -# include -# else +# include +# else # include -# endif +# endif #endif // First, define library-dependent types (complex, FFT plan) @@ -48,9 +48,9 @@ namespace warpx::utils # else using Complex = fftw_complex; # endif -# else +# else using Complex = std::complex; -# endif +# endif #endif } -#endif // WARPX_UTILS_MATH_COMPLEX__H_ \ No newline at end of file +#endif // WARPX_UTILS_MATH_COMPLEX__H_ diff --git a/Source/WarpX.H b/Source/WarpX.H index 6c15d4f63a2..3c7d7fe45ef 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -335,7 +335,7 @@ public: static bool safe_guard_cells; - //enable radiations in the evolve + //enable radiations in the evolve bool m_do_radiation_flag; //! With mesh refinement, particles located inside a refinement patch, but within diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 130e51c970e..55bf54a2b11 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1919,7 +1919,7 @@ WarpX::BackwardCompatibility () "collisions.ncollisions is ignored. Just use particles.collision_names please.", ablastr::warn_manager::WarnPriority::low); } - + const ParmParse pp_lasers("lasers"); int nlasers; if (pp_lasers.query("nlasers", nlasers)){ From be0ae475c1e351caff84af0173af395a9284e137 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 15:07:27 +0100 Subject: [PATCH 14/60] use false insteand of 0 for bool --- Source/WarpX.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 55bf54a2b11..d067fdb99c1 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -215,7 +215,7 @@ bool WarpX::do_device_synchronize = true; bool WarpX::do_device_synchronize = false; #endif -bool m_do_radiation_flag = 0; +bool m_do_radiation_flag = false; WarpX* WarpX::m_instance = nullptr; From 2e4b9f8ddcc59067979fe4c15e25c091f7b3296c Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 15:10:29 +0100 Subject: [PATCH 15/60] remove radiation flag from warpx class --- Source/WarpX.H | 4 ---- Source/WarpX.cpp | 2 -- 2 files changed, 6 deletions(-) diff --git a/Source/WarpX.H b/Source/WarpX.H index 3c7d7fe45ef..c644aa67982 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -334,10 +334,6 @@ public: static bool do_device_synchronize; static bool safe_guard_cells; - - //enable radiations in the evolve - bool m_do_radiation_flag; - //! With mesh refinement, particles located inside a refinement patch, but within //! #n_field_gather_buffer cells of the edge of the patch, will gather the fields //! from the lower refinement level instead of the refinement patch itself diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index d067fdb99c1..e862391383e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -215,8 +215,6 @@ bool WarpX::do_device_synchronize = true; bool WarpX::do_device_synchronize = false; #endif -bool m_do_radiation_flag = false; - WarpX* WarpX::m_instance = nullptr; void WarpX::MakeWarpX () From 091b68bc8ca9d3ada123cf00819954087507acf6 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 15:14:44 +0100 Subject: [PATCH 16/60] remove Utils/Math/Complex.H --- Source/Particles/Radiation/RadiationHandler.H | 1 - Source/Utils/Math/Complex.H | 56 ------------------- 2 files changed, 57 deletions(-) delete mode 100644 Source/Utils/Math/Complex.H diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 4d195005879..469b908cc67 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -8,7 +8,6 @@ #ifndef WARPX_PARTICLES_RADIATION_H_ #define WARPX_PARTICLES_RADIATION_H_ -#include "Utils/Math/Complex.H" #include "Particles/Pusher/GetAndSetPosition.H" diff --git a/Source/Utils/Math/Complex.H b/Source/Utils/Math/Complex.H deleted file mode 100644 index 651dc3449dc..00000000000 --- a/Source/Utils/Math/Complex.H +++ /dev/null @@ -1,56 +0,0 @@ -/* Copyright 2022 Thomas Clark - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ - -#ifndef WARPX_UTILS_MATH_COMPLEX__H_ -#define WARPX_UTILS_MATH_COMPLEX__H_ - -#include -#include - -#if defined(AMREX_USE_CUDA) -# include -#elif defined(AMREX_USE_HIP) -# include - -#else -# if WARPX_USE_PSATD -# include -# else -# include - -# endif -#endif - - // First, define library-dependent types (complex, FFT plan) -namespace warpx::utils -{ - /** Complex type for FFT, depends on FFT library */ -#if defined(AMREX_USE_CUDA) -# ifdef AMREX_USE_FLOAT - using Complex = cuComplex; -# else - using Complex = cuDoubleComplex; -# endif -#elif defined(AMREX_USE_HIP) -# ifdef AMREX_USE_FLOAT - using Complex = float2; -# else - using Complex = double2; -# endif -#else -#if WARPX_USE_PSATD -# ifdef AMREX_USE_FLOAT - using Complex = fftwf_complex; -# else - using Complex = fftw_complex; -# endif -# else - using Complex = std::complex; -# endif -#endif -} -#endif // WARPX_UTILS_MATH_COMPLEX__H_ From 4c4397f604bc321b320377c75cc8ee882ab11bb7 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 15:31:16 +0100 Subject: [PATCH 17/60] move WarpX_Complex.H into ablastr --- .../PsatdAlgorithmComoving.cpp | 6 ++-- .../PsatdAlgorithmFirstOrder.cpp | 4 ++- .../PsatdAlgorithmJConstantInTime.cpp | 3 +- .../PsatdAlgorithmJLinearInTime.cpp | 4 ++- .../SpectralAlgorithms/PsatdAlgorithmPml.cpp | 4 ++- .../SpectralBaseAlgorithm.H | 5 +-- .../SpectralBaseAlgorithm.cpp | 4 ++- .../SpectralSolver/SpectralFieldData.H | 4 +-- .../SpectralSolver/SpectralKSpace.H | 4 +-- .../LaserProfileFieldFunction.cpp | 1 - .../LaserProfileFromFile.cpp | 3 +- .../LaserProfileGaussian.cpp | 5 +-- .../Particles/Deposition/ChargeDeposition.H | 5 ++- .../Particles/Deposition/CurrentDeposition.H | 9 +++-- .../Deposition/SharedDepositionUtils.H | 4 ++- Source/Particles/Gather/FieldGather.H | 4 ++- Source/Utils/WarpX_Complex.H | 32 ----------------- Source/ablastr/math/Complex.H | 36 +++++++++++++++++++ 18 files changed, 82 insertions(+), 55 deletions(-) delete mode 100644 Source/Utils/WarpX_Complex.H create mode 100644 Source/ablastr/math/Complex.H 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..3e564f74b31 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 ablsatr::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/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/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/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index d0db678dfda..bbc349be243 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 @@ -51,6 +52,7 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, long load_balance_costs_update_algo) { using namespace amrex; + using namespace ablastr::math; #if !defined(AMREX_USE_GPU) amrex::ignore_unused(cost, load_balance_costs_update_algo); @@ -260,6 +262,7 @@ void doChargeDepositionSharedShapeN (const GetParticlePosition& GetPositio ) { using namespace amrex; + using namespace ablastr::math; const auto *permutation = a_bins.permutationPtr(); diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index ac377a6799d..94cfb0d3d7a 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,7 @@ void doDepositionShapeN (const GetParticlePosition& GetPosition, long load_balance_costs_update_algo) { using namespace amrex::literals; + using namespace ablastr::math; #if !defined(WARPX_DIM_RZ) amrex::ignore_unused(n_rz_azimuthal_modes); @@ -599,6 +601,7 @@ void doEsirkepovDepositionShapeN (const GetParticlePosition& GetPosition, { using namespace amrex; using namespace amrex::literals; + using namespace ablastr::math; #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..f3da11e94a9 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 @@ -83,6 +84,7 @@ void depositComponent (const GetParticlePosition& GetPosition, const int zdir, const int NODE, const int CELL, const int dir) { using namespace amrex::literals; + using namespace ablastr::math; #if !defined(WARPX_DIM_RZ) amrex::ignore_unused(n_rz_azimuthal_modes); 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/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..e6cf1db1ee6 --- /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_ From 5ca02819b2cfc4528f89695e7d683e326bd99682 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 15:59:41 +0100 Subject: [PATCH 18/60] fix bugs --- Source/Laser/LaserProfiles.H | 6 +++--- Source/Particles/Deposition/ChargeDeposition.H | 5 +++++ Source/Particles/Deposition/CurrentDeposition.H | 4 ++++ Source/Particles/Radiation/RadiationHandler.cpp | 17 +++++++++-------- 4 files changed, 21 insertions(+), 11 deletions(-) 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/Particles/Deposition/ChargeDeposition.H b/Source/Particles/Deposition/ChargeDeposition.H index bbc349be243..128a0edabc4 100644 --- a/Source/Particles/Deposition/ChargeDeposition.H +++ b/Source/Particles/Deposition/ChargeDeposition.H @@ -52,7 +52,10 @@ void doChargeDepositionShapeN (const GetParticlePosition& GetPosition, long load_balance_costs_update_algo) { 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); @@ -262,7 +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 94cfb0d3d7a..c42708fcf0d 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -72,7 +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); @@ -601,7 +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/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 02d58e0e40e..d8d7b5d484e 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -12,10 +12,11 @@ #include #include - +#include #include -using namespace warpx::utils; + +using namespace ablastr::math; RadiationHandler::RadiationHandler() { @@ -164,12 +165,12 @@ void RadiationHandler::add_radiation_contribution Complex Term_z= q*dt/(16*pow(ablastr::constant::math::pi,3)*ablastr::constant::SI::ep0*ablastr::constant::SI::c)*ncrossncrossBetapointz/std::pow(un_betan,2)*eiomega; //Add the contributions - m_radiation_data[i_x,i_y,i_om,0]+=std::real(Term_x); - m_radiation_data[i_x,i_y,i_om,1]+=std::imag(Term_x); - m_radiation_data[i_x,i_y,i_om,2]+=std::real(Term_y); - m_radiation_data[i_x,i_y,i_om,3]+=std::imag(Term_y); - m_radiation_data[i_x,i_y,i_om,4]+=std::real(Term_z); - m_radiation_data[i_x,i_y,i_om,5]+=std::imag(Term_z); + m_radiation_data[i_x,i_y,i_om,0]+=Term_x.m_real; + m_radiation_data[i_x,i_y,i_om,1]+=Term_x.m_imag; + m_radiation_data[i_x,i_y,i_om,2]+=Term_y.m_real; + m_radiation_data[i_x,i_y,i_om,3]+=Term_y.m_imag; + m_radiation_data[i_x,i_y,i_om,4]+=Term_z.m_real; + m_radiation_data[i_x,i_y,i_om,5]+=Term_z.m_imag; } } From 05eec9955eab3179aa4ba46f2014eec81bbf87d6 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 16:11:03 +0100 Subject: [PATCH 19/60] remove unused function --- Source/Particles/Radiation/RadiationHandler.cpp | 1 - Source/Particles/Sorting/SortingUtils.H | 8 -------- 2 files changed, 9 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index d8d7b5d484e..970f85776c0 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -233,7 +233,6 @@ void RadiationHandler::add_detector det_bornes[idim][0]=center[idim]-std::tan(m_theta_range[idim]/2)*m_det_distance; det_bornes[idim][1]=center[idim]+std::tan(m_theta_range[idim]/2)*m_det_distance; } - //fillWithConsecutiveReal(pos_det_x) //pos_det_y //omega_calc diff --git a/Source/Particles/Sorting/SortingUtils.H b/Source/Particles/Sorting/SortingUtils.H index 076c6a21e0c..1bb398256cf 100644 --- a/Source/Particles/Sorting/SortingUtils.H +++ b/Source/Particles/Sorting/SortingUtils.H @@ -21,14 +21,6 @@ */ void fillWithConsecutiveIntegers( amrex::Gpu::DeviceVector& v ); -/** \brief Fill the elements of the input vector with consecutive integer, - * starting from 0 - * - * \param[inout] v Vector of reals, to be filled by this routine - */ - -void fillWithConsecutiveReal( amrex::Gpu::DeviceVector& v ); - /** \brief Find the indices that would reorder the elements of `predicate` * so that the elements with non-zero value precede the other elements From 8ca7b6a31ec2071c1b466794d09d3c56e16ecc24 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 16:12:44 +0100 Subject: [PATCH 20/60] removed comments --- Source/Particles/MultiParticleContainer.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index aaf3bcaf315..1194cf86b29 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -122,13 +122,11 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) // Setup particle collisions collisionhandler = std::make_unique(this); - amrex::Print() << "wEEEEEEEESSSSSSSHHHHH" << allcontainers.size() << "\n"; //Initialization of the radiation for (auto& s : allcontainers) { amrex::Print() << s->has_radiation() <<"jzfnfhjzbefhezbfhzebfzehb" << "\n"; if (s->has_radiation()) { m_at_least_one_has_radiation = true; - amrex::Print() << "Hey ! I'm here" << "\n"; m_p_radiation_handler = std::make_unique(); break; } From 706ca3be953766a6307621450b56611ec0d0f979 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 16:23:53 +0100 Subject: [PATCH 21/60] using forward declaration for radiation module --- Source/Particles/MultiParticleContainer.H | 4 ++-- Source/Particles/MultiParticleContainer.cpp | 3 +++ Source/Particles/Radiation/RadiationHandler.H | 20 +++++++++---------- .../Radiation/RadiationHandler_fwd.H | 17 ++++++++++++++++ 4 files changed, 32 insertions(+), 12 deletions(-) create mode 100644 Source/Particles/Radiation/RadiationHandler_fwd.H diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 016b4e84380..ecb6beff344 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -20,7 +20,7 @@ # include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper_fwd.H" #endif #include "PhysicalParticleContainer.H" -#include "Particles/Radiation/RadiationHandler.H" +#include "Particles/Radiation/RadiationHandler_fwd.H" #include "Utils/TextMsg.H" #include "Utils/WarpXConst.H" #include "WarpXParticleContainer.H" @@ -69,7 +69,7 @@ public: MultiParticleContainer (amrex::AmrCore* amr_core); - ~MultiParticleContainer() = default; + ~MultiParticleContainer(); MultiParticleContainer (MultiParticleContainer const &) = delete; MultiParticleContainer& operator= (MultiParticleContainer const & ) = delete; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 1194cf86b29..66ec7e118f0 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" @@ -133,6 +134,8 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) } } +MultiParticleContainer::~MultiParticleContainer() {} + void MultiParticleContainer::ReadParameters () { diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 469b908cc67..016149e5d9f 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -6,24 +6,24 @@ * License: BSD-3-Clause-LBNL */ -#ifndef WARPX_PARTICLES_RADIATION_H_ -#define WARPX_PARTICLES_RADIATION_H_ -#include "Particles/Pusher/GetAndSetPosition.H" +#ifndef WARPX_PARTICLES_RADIATION_H +#define WARPX_PARTICLES_RADIATION_H + +#include "RadiationHandler_fwd.H" +#include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/Sorting/SortingUtils.H" +#include "WarpX.H" #include #include -#include -#include -#include #include -#include "WarpX.H" -#include "Particles/Sorting/SortingUtils.H" + #include #include -/* \brief CollisionHandler is a light weight class that contains the +/* \brief CollisionHandler is a class that contains the * calculation of radiation for particles at each frequencies and angles. */ class RadiationHandler @@ -82,4 +82,4 @@ private: amrex::Vector omega_calc; }; -#endif // WARPX_PARTICLES_RADIATION_H_ +#endif // WARPX_PARTICLES_RADIATION_H diff --git a/Source/Particles/Radiation/RadiationHandler_fwd.H b/Source/Particles/Radiation/RadiationHandler_fwd.H new file mode 100644 index 00000000000..6b27274272b --- /dev/null +++ b/Source/Particles/Radiation/RadiationHandler_fwd.H @@ -0,0 +1,17 @@ + +/* Copyright 2023 Thomas Clark, Luca Fedeli + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_PARTICLES_RADIATION_FWD_H +#define WARPX_PARTICLES_RADIATION_FWD_H + +/* \brief CollisionHandler is a class that contains the + * calculation of radiation for particles at each frequencies and angles. + */ +class RadiationHandler; + +#endif // WARPX_PARTICLES_RADIATION_FWD_H From dc2641240ec3271309bfe0ef580a7dcf159a0be9 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 16:26:37 +0100 Subject: [PATCH 22/60] remove unused includes --- Source/Particles/Radiation/RadiationHandler.H | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 016149e5d9f..44b27bd8064 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -13,12 +13,10 @@ #include "Particles/Pusher/GetAndSetPosition.H" #include "Particles/Sorting/SortingUtils.H" -#include "WarpX.H" +#include #include #include -#include - #include #include From 0c2f68e1b12d137f7e4de368bc0e59938e25fbbe Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 16:39:06 +0100 Subject: [PATCH 23/60] remove text message --- Source/Particles/MultiParticleContainer.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 66ec7e118f0..3ae58e4e8a1 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -125,7 +125,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) collisionhandler = std::make_unique(this); //Initialization of the radiation for (auto& s : allcontainers) { - amrex::Print() << s->has_radiation() <<"jzfnfhjzbefhezbfhzebfzehb" << "\n"; if (s->has_radiation()) { m_at_least_one_has_radiation = true; m_p_radiation_handler = std::make_unique(); From 6e50dc4798f5056040c6e485b5c731be00ef6ae8 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 15 Dec 2023 17:37:33 +0100 Subject: [PATCH 24/60] refactoring --- Source/Evolve/WarpXEvolve.cpp | 4 +- Source/Particles/MultiParticleContainer.H | 3 +- Source/Particles/MultiParticleContainer.cpp | 48 ++++--------------- .../Particles/PhysicalParticleContainer.cpp | 12 ++--- Source/Particles/Radiation/RadiationHandler.H | 4 +- .../Particles/Radiation/RadiationHandler.cpp | 29 +++++------ Source/Particles/WarpXParticleContainer.H | 13 +++++ Source/Particles/WarpXParticleContainer.cpp | 41 ++++++++++++++++ 8 files changed, 85 insertions(+), 69 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index 2d5101d7cdc..f7cef17c8e2 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -428,8 +428,8 @@ WarpX::OneStep_nosub (Real cur_time) ExecutePythonCallback("particlescraper"); ExecutePythonCallback("beforedeposition"); - //Save particle old momentum in a attribute - mypc->keepoldmomentum(); + //Save particles' old momenta + mypc->RecordOldMomenta(); PushParticlesandDeposit(cur_time); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index ecb6beff344..85b95894513 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -332,8 +332,7 @@ public: void doRadiation(amrex::Real dt, amrex::Real cur_time); - void keepoldmomentum(); - void keepoldmomentum_p(std::unique_ptr& pc); + void RecordOldMomenta(); void Dump_radiations(); [[nodiscard]] int getSpeciesID (std::string product_str) const; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 3ae58e4e8a1..748446525be 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -127,7 +127,13 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) for (auto& s : allcontainers) { if (s->has_radiation()) { m_at_least_one_has_radiation = true; - m_p_radiation_handler = std::make_unique(); + + const auto& lev_0_geom = amr_core->Geom(0); + auto center = amrex::Array{}; + for(int idim=0; idim(center); break; } } @@ -923,46 +929,13 @@ MultiParticleContainer::doFieldIonization (int lev, } } -void MultiParticleContainer::keepoldmomentum(){ +void MultiParticleContainer::RecordOldMomenta(){ for (auto& pc : allcontainers) { if (pc->has_radiation()){ - keepoldmomentum_p(pc); - } - } -}void MultiParticleContainer::keepoldmomentum_p - (std::unique_ptr& pc){ - const auto level0=0; - for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { - auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); - auto& part=pc->GetParticles(level0)[index]; - long const np = pti.numParticles(); - auto& soa = part.GetStructOfArrays(); - - //Load the momentums - amrex::ParticleReal* ux = soa.GetRealData(PIdx::ux).data(); - amrex::ParticleReal* uy = soa.GetRealData(PIdx::uy).data(); - amrex::ParticleReal* uz = soa.GetRealData(PIdx::uz).data(); - - //Finding the good attribute index - int index_name_x=pc->GetRealCompIndex("prev_u_x"); - int index_name_y=pc->GetRealCompIndex("prev_u_y"); - int index_name_z=pc->GetRealCompIndex("prev_u_z"); - - auto* p_ux = soa.GetRealData(index_name_x).data(); - auto* p_uy = soa.GetRealData(index_name_y).data(); - auto* p_uz = soa.GetRealData(index_name_z).data(); - - amrex::ParallelFor(np, - [=] AMREX_GPU_DEVICE(int ip) - { - //Putting them in an attribute - p_ux[ip] = ux[ip]; - p_uy[ip] = uy[ip]; - p_uz[ip] = uz[ip]; - - }); + pc->RecordOldMomenta(); } } +} void MultiParticleContainer::doCollisions ( Real cur_time, amrex::Real dt ) @@ -1009,7 +982,6 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt, amrex::Real cur_ if (m_at_least_one_has_radiation){ for (auto& pc : allcontainers) { if (pc->has_radiation()){ - m_p_radiation_handler->add_radiation_contribution(dt,pc,cur_time); } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index e386bdbbe46..645f90f0937 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -347,15 +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 previous momentum attribute + + //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){ - amrex::Print() << "i'm here"; - //enable the radiations in the evolve loop - AddRealComp("prev_u_x"); - AddRealComp("prev_u_y"); - AddRealComp("prev_u_z"); - } + 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/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 44b27bd8064..37f1cddf2eb 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -27,17 +27,17 @@ class RadiationHandler { public: - RadiationHandler (); + 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 add_detector(); void gather_and_write_radiation(const std::string& filename); void Integral_overtime(); private: + void add_detector(const amrex::Array& center); // Frequency range amrex::Vector m_omega_range; diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 970f85776c0..ae600a8c7a0 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -18,7 +18,7 @@ using namespace ablastr::math; -RadiationHandler::RadiationHandler() +RadiationHandler::RadiationHandler(const amrex::Array& center) { // Read in radiation input const amrex::ParmParse pp_radiations("radiations"); @@ -43,7 +43,8 @@ RadiationHandler::RadiationHandler() pp_radiations.getarr("detector_number_points", m_det_pts,0,2); pp_radiations.getarr("detector_direction", m_det_direction,0,3); pp_radiations.query("detector_distance", m_det_distance); - add_detector(); + + add_detector(center); /* Initialize the Fab with the field in function of angle and frequency */ //int numcomps = 4; //BaseFab fab(bx,numcomps); @@ -93,13 +94,13 @@ void RadiationHandler::add_radiation_contribution auto& part=pc->GetParticles(level0)[index]; auto& soa = part.GetStructOfArrays(); - const auto prev_u_x_idx =pc->GetRealCompIndex("prev_u_x"); - const auto prev_u_y_idx =pc->GetRealCompIndex("prev_u_y"); - const auto prev_u_z_idx =pc->GetRealCompIndex("prev_u_z"); + const auto old_u_x_idx =pc->GetRealCompIndex("old_u_x"); + const auto old_u_y_idx =pc->GetRealCompIndex("old_u_y"); + const auto old_u_z_idx =pc->GetRealCompIndex("old_u_z"); - auto* p_ux_old = soa.GetRealData(prev_u_x_idx).data(); - auto* p_uy_old = soa.GetRealData(prev_u_y_idx).data(); - auto* p_uz_old = soa.GetRealData(prev_u_z_idx).data(); + auto* p_ux_old = soa.GetRealData(old_u_x_idx).data(); + auto* p_uy_old = soa.GetRealData(old_u_y_idx).data(); + auto* p_uz_old = soa.GetRealData(old_u_z_idx).data(); auto GetPosition = GetParticlePosition(pti); amrex::ParticleReal const q = pc->getCharge(); @@ -200,10 +201,8 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) of.close(); } } -void RadiationHandler::add_detector - (){ - const auto level0=0; - amrex::Geometry const& geom = WarpX::GetInstance().Geom(level0); +void RadiationHandler::add_detector(const amrex::Array& center){ + //Calculation of angle resolution m_d_theta.resize(2); for(int i=0; i<2; i++){ @@ -211,11 +210,7 @@ void RadiationHandler::add_detector } //Set the resolution of the detector m_d_d.resize(3); - amrex::Vector center; - center.resize(3); - for(int idim=0; idim Date: Fri, 15 Dec 2023 17:42:39 +0100 Subject: [PATCH 25/60] fixed bug --- Source/ablastr/math/Complex.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/ablastr/math/Complex.H b/Source/ablastr/math/Complex.H index e6cf1db1ee6..705e5182bbb 100644 --- a/Source/ablastr/math/Complex.H +++ b/Source/ablastr/math/Complex.H @@ -31,6 +31,6 @@ namespace ablastr::math static_assert(sizeof(Complex) == sizeof(amrex::Real[2]), "Unexpected complex type."); -}; +} #endif //ABLASTR_MATH_COMPLEX_H_ From 2282ebda111ff12b8eb00c8d9c2099c20a0db99e Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Mon, 18 Dec 2023 16:04:51 +0100 Subject: [PATCH 26/60] cleaning --- Source/Particles/Radiation/RadiationHandler.H | 39 +---- .../Particles/Radiation/RadiationHandler.cpp | 164 +++++++++--------- Source/ablastr/math/Utils.H | 30 ++++ 3 files changed, 123 insertions(+), 110 deletions(-) create mode 100644 Source/ablastr/math/Utils.H diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 37f1cddf2eb..6dd75d76fed 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -11,9 +11,6 @@ #include "RadiationHandler_fwd.H" -#include "Particles/Pusher/GetAndSetPosition.H" -#include "Particles/Sorting/SortingUtils.H" - #include #include #include @@ -32,50 +29,28 @@ public: /* 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(); + void Integral_overtime(); private: - void add_detector(const amrex::Array& center); // Frequency range - amrex::Vector m_omega_range; + amrex::Array m_omega_range; // Dimensions of the detector - amrex::Vector m_theta_range; + amrex::Array m_theta_range; - int m_omega_points; - - //put a detector - int m_get_a_detector; - // Resolution of the detector - amrex::Vector m_d_theta; - double m_d_omega; - //The vector with the resolution of the detector - amrex::Vector m_d_det; amrex::Real m_det_distance; + amrex::Array m_det_pts; + amrex::Array m_det_direction; + amrex::Array m_det_orientation; - //Define the Fab with the datas - int ncomp; - - amrex::Vector m_det_pts; - amrex::Vector m_det_direction; - amrex::Vector m_pos_det; - amrex::Vector n; - - //resolution of the grid of the detectir - amrex::Vector m_d_d; - - // amrex::Gpu::DeviceVector m_radiation_data; - amrex::Gpu::DeviceVector m_radiation_calculation; - // amrex::Real dephas; - // - amrex::Vector> det_bornes; amrex::Vector> det_pos; amrex::Vector omega_calc; diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index ae600a8c7a0..067dc66d0bf 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -7,72 +7,112 @@ #include "RadiationHandler.H" +#include "Particles/Pusher/GetAndSetPosition.H" + #include "Utils/TextMsg.H" #include #include #include +#include #include using namespace ablastr::math; -RadiationHandler::RadiationHandler(const amrex::Array& center) +namespace { - // Read in radiation input - const amrex::ParmParse pp_radiations("radiations"); - - // Verify if there is a detector - pp_radiations.query("put_a_detector", m_get_a_detector); + auto compute_detector_positions( + const amrex::Array& center, + const amrex::Array& direction, + amrex::Real distance, + const amrex::Array& orientation, + amrex::Array& det_points, + amrex::Array& theta_range + ) + { + + } +} +RadiationHandler::RadiationHandler(const amrex::Array& center) +{ #if defined WARPX_DIM_RZ WARPX_ABORT_WITH_MESSAGE("Radiation is not supported yet with RZ."); #endif - //Compute the radiation - if(m_get_a_detector){ - //Resolution in frequency of the detector - pp_radiations.queryarr("omega_range", m_omega_range); - pp_radiations.query("omega_points", m_omega_points); - //Angle theta AND phi - pp_radiations.queryarr("angle_aperture", m_theta_range); - - //Type of detector - pp_radiations.getarr("detector_number_points", m_det_pts,0,2); - pp_radiations.getarr("detector_direction", m_det_direction,0,3); - pp_radiations.query("detector_distance", m_det_distance); - - add_detector(center); - /* Initialize the Fab with the field in function of angle and frequency */ - //int numcomps = 4; - //BaseFab fab(bx,numcomps); - // Box of angle and frequency - //const amrex::Box detect_box({0,0,0}, {m_det_pts[0], m_det_pts[1], m_omega_points}); - ncomp = 6; //Real and complex part - //amrex::FArrayBox fab_detect(detect_box, ncomp); - //fab_detect.setval(0,0) - m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); - m_radiation_calculation = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points); - - det_pos.resize(2); - det_pos[0].resize(m_det_pts[0]); - det_pos[1].resize(m_det_pts[1]); - omega_calc.resize(m_omega_points); - - //Initialization : - for(int i=0; i(2); + pp_radiation.getarr("omega_range", omega_range); + std::copy(omega_range.begin(), omega_range.end(), m_omega_range.begin()); + pp_radiation.get("omega_points", m_omega_points); + + //Angle theta AND phi + auto theta_range std::vector(2); + pp_radiation.get("angle_aperture", theta_range); + std::copy(theta_range.begin(), theta_range.end(), m_theta_range.begin()); + + //Detector parameters + auto det_pts std::vector(2); + pp_radiation.getarr("detector_number_points", det_pts); + std::copy(det_pts.begin(), det_pts.end(), m_det_pts.begin()); + + auto det_direction std::vector(3); + pp_radiation.getarr("detector_direction", det_direction); + std::copy(det_direction.begin(), det_direction.end(), m_det_direction.begin()); + + auto det_orientation std::vector(3); + pp_radiation.getarr("detector_orientation", det_orientation); + std::copy(det_orientation.begin(), det_orientation.end(), m_det_orientation.begin()); + + pp_radiation.get("detector_distance", m_det_distance); + + //Calculation of angle resolution + const auto d_theta = amrex::Array{ + 2*m_theta_range[0]/m_det_pts[0], + 2*m_theta_range[1]/m_det_pts[1]}; + + //Set the resolution of the detector + m_d_d.resize(3); + + for(int idi = 0; idi<3; ++idi){ + if(m_det_direction[idi]==1){ + m_d_d[(idi+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); + m_d_d[(idi+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); + } + + m_d_omega=m_omega_range[1]-m_omega_range[0]/static_cast(m_omega_points); + + //Calculate the sides of the detector + det_bornes.resize(2); + det_bornes[0].resize(2); + det_bornes[1].resize(2); + + for(int idim = 0; idim<2; ++idim){ + det_bornes[idim][0]=center[idim]-std::tan(m_theta_range[idim]/2)*m_det_distance; + det_bornes[idim][1]=center[idim]+std::tan(m_theta_range[idim]/2)*m_det_distance; + } + //pos_det_y + //omega_calc } + constexpr auto ncomp = 6; + m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); + + det_pos.resize(2); + det_pos[0].resize(m_det_pts[0]); + det_pos[1].resize(m_det_pts[1]); + omega_calc.resize(m_omega_points); + + ablastr::math::LinSpaceFill(det_bornes[0][0], det_bornes[0][1], m_det_pts[0], det_pos[0].being()); + ablastr::math::LinSpaceFill(det_bornes[1][0], det_bornes[1][1], m_det_pts[1], det_pos[1].being()); + + ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, omega_calc.being()); } void RadiationHandler::add_radiation_contribution @@ -201,38 +241,6 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) of.close(); } } -void RadiationHandler::add_detector(const amrex::Array& center){ - - //Calculation of angle resolution - m_d_theta.resize(2); - for(int i=0; i<2; i++){ - m_d_theta[i] = 2*m_theta_range[i]/static_cast(m_det_pts[i]); - } - //Set the resolution of the detector - m_d_d.resize(3); - - for(int idi = 0; idi<3; ++idi){ - if(m_det_direction[idi]==1){ - m_d_d[(idi+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); - m_d_d[(idi+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); - } - - m_d_omega=m_omega_range[1]-m_omega_range[0]/static_cast(m_omega_points); - - //Calculate the sides of the detector - det_bornes.resize(2); - det_bornes[0].resize(2); - det_bornes[1].resize(2); - - for(int idim = 0; idim<2; ++idim){ - det_bornes[idim][0]=center[idim]-std::tan(m_theta_range[idim]/2)*m_det_distance; - det_bornes[idim][1]=center[idim]+std::tan(m_theta_range[idim]/2)*m_det_distance; - } - //pos_det_y - //omega_calc - - } -} void RadiationHandler::Integral_overtime() { diff --git a/Source/ablastr/math/Utils.H b/Source/ablastr/math/Utils.H new file mode 100644 index 00000000000..a94e1da484d --- /dev/null +++ b/Source/ablastr/math/Utils.H @@ -0,0 +1,30 @@ +/* Copyright 2023 Thomas Clark, Luca Fedeli + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ +#ifndef ABLASTR_MATH_UTILS_H_ +#define ABLASTR_MATH_UTILS_H_ + +#include +#include + +namespace ablastr::math +{ + template + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + void LinSpaceFill(T start, T end, int size, Iter& it) + { + const auto delta = static_cast((end - start)/(size-1)); + + for (int i = 0; i < size-1; ++i) + { + *(it++) = start + static_cast(i*delta); + } + *it = end; + } + +} + +#endif //ABLASTR_MATH_UTILS_H_ From 99c42abfa3ac053bf486e247a870898946af4017 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 19 Dec 2023 13:57:25 +0100 Subject: [PATCH 27/60] continue cleaning --- Source/Particles/Radiation/RadiationHandler.H | 14 +- .../Particles/Radiation/RadiationHandler.cpp | 175 ++++++++++-------- 2 files changed, 106 insertions(+), 83 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 6dd75d76fed..316f08c4f54 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -11,6 +11,9 @@ #include "RadiationHandler_fwd.H" +#include + + #include #include #include @@ -38,6 +41,8 @@ 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; @@ -47,12 +52,11 @@ private: amrex::Array m_det_direction; amrex::Array m_det_orientation; - amrex::Gpu::DeviceVector m_radiation_data; - - amrex::Real dephas; + amrex::Gpu::DeviceVector det_pos_x; + amrex::Gpu::DeviceVector det_pos_y; + amrex::Gpu::DeviceVector det_pos_z; - amrex::Vector> det_pos; - amrex::Vector omega_calc; + amrex::Gpu::DeviceVector m_radiation_data; }; #endif // WARPX_PARTICLES_RADIATION_H diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 067dc66d0bf..35b777e28f3 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -8,15 +8,13 @@ #include "RadiationHandler.H" #include "Particles/Pusher/GetAndSetPosition.H" - #include "Utils/TextMsg.H" -#include - #include -#include #include +#include + #include using namespace ablastr::math; @@ -26,13 +24,70 @@ namespace auto compute_detector_positions( const amrex::Array& center, const amrex::Array& direction, - amrex::Real distance, + const amrex::Real distance, const amrex::Array& orientation, - amrex::Array& det_points, - amrex::Array& theta_range - ) + const amrex::Array& det_points, + const amrex::Array& theta_range) { - + 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); + ablastr::math::LinSpaceFill(-ul,ul, det_points[0], us.being()); + + auto vs = amrex::Vector(det_points[1]); + const auto vlim = distance*std::tan(theta_range[1]*0.5_rt); + ablastr::math::LinSpaceFill(-vlim, vlim, det_points[1], vs.being()); + + const auto how_many = det_points[0]*det_points[1]; + + auto det_x = amrex::GPU::DeviceVector(how_many); + auto det_y = amrex::GPU::DeviceVector(how_many); + auto det_z = amrex::GPU::DeviceVector(how_many); + + 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) + { + det_x[i*det_points[1] + j] = center[0] + distance * norm_direction[0] + us[i]*u[0] + vs[j]*v[0]; + det_x[i*det_points[1] + j] = center[1] + distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; + det_x[i*det_points[1] + j] = center[2] + distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; + } + } + + amrex::GPU::synchronize(); + + return {det_x, det_y, det_z} + } } @@ -72,57 +127,30 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) pp_radiation.get("detector_distance", m_det_distance); - //Calculation of angle resolution - const auto d_theta = amrex::Array{ - 2*m_theta_range[0]/m_det_pts[0], - 2*m_theta_range[1]/m_det_pts[1]}; - - //Set the resolution of the detector - m_d_d.resize(3); - - for(int idi = 0; idi<3; ++idi){ - if(m_det_direction[idi]==1){ - m_d_d[(idi+1)%3]=2*m_det_distance*tan(m_d_theta[0]/2); - m_d_d[(idi+2)%3]=2*m_det_distance*tan(m_d_theta[1]/2); - } - - m_d_omega=m_omega_range[1]-m_omega_range[0]/static_cast(m_omega_points); - - //Calculate the sides of the detector - det_bornes.resize(2); - det_bornes[0].resize(2); - det_bornes[1].resize(2); - - for(int idim = 0; idim<2; ++idim){ - det_bornes[idim][0]=center[idim]-std::tan(m_theta_range[idim]/2)*m_det_distance; - det_bornes[idim][1]=center[idim]+std::tan(m_theta_range[idim]/2)*m_det_distance; - } - //pos_det_y - //omega_calc - - } - constexpr auto ncomp = 6; - m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); + [det_pos_x, det_pos_y, det_pos_z] = compute_detector_positions( + center, det_direction, m_det_distance, + m_det_orientation, m_det_pts, m_theta_range); - det_pos.resize(2); - det_pos[0].resize(m_det_pts[0]); - det_pos[1].resize(m_det_pts[1]); - omega_calc.resize(m_omega_points); - - ablastr::math::LinSpaceFill(det_bornes[0][0], det_bornes[0][1], m_det_pts[0], det_pos[0].being()); - ablastr::math::LinSpaceFill(det_bornes[1][0], det_bornes[1][1], m_det_pts[1], det_pos[1].being()); + constexpr auto ncomp = 3; + m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); + m_omegas = amrex::Gpu::DeviceVector(m_omega_points); ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, omega_calc.being()); } + void RadiationHandler::add_radiation_contribution - (const amrex::Real dt, std::unique_ptr& pc, amrex::Real current_time){ - const auto level0=0; + (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 + #pragma omp parallel #endif -{ - for (WarpXParIter pti(*pc, level0); pti.isValid(); ++pti) { + { + for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) + { long const np = pti.numParticles(); auto& attribs = pti.GetAttribs(); auto& p_w = attribs[PIdx::w]; @@ -131,39 +159,23 @@ void RadiationHandler::add_radiation_contribution auto& p_uz = attribs[PIdx::uz]; auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); - auto& part=pc->GetParticles(level0)[index]; + auto& part=pc->GetParticles(lev)[index]; auto& soa = part.GetStructOfArrays(); - const auto old_u_x_idx =pc->GetRealCompIndex("old_u_x"); - const auto old_u_y_idx =pc->GetRealCompIndex("old_u_y"); - const auto old_u_z_idx =pc->GetRealCompIndex("old_u_z"); - - auto* p_ux_old = soa.GetRealData(old_u_x_idx).data(); - auto* p_uy_old = soa.GetRealData(old_u_y_idx).data(); - auto* p_uz_old = soa.GetRealData(old_u_z_idx).data(); + auto* p_ux_old = soa.GetRealData(pc->GetRealCompIndex("old_u_x")).data(); + auto* p_uy_old = soa.GetRealData(pc->GetRealCompIndex("old_u_y")).data(); + auto* p_uz_old = soa.GetRealData(pc->GetRealCompIndex("old_u_z")).data(); auto GetPosition = GetParticlePosition(pti); amrex::ParticleReal const q = pc->getCharge(); - amrex::ParallelFor(np, - [=] AMREX_GPU_DEVICE(int ip) - { amrex::ParticleReal xp, yp, zp; - - GetPosition.AsStored(ip,xp, yp, zp); - amrex::GpuArray Part_pos{xp,yp,zp}; - - amrex::ParticleReal p_ux_old_unit=p_ux_old[ip]; - amrex::ParticleReal p_uy_old_unit=p_uy_old[ip]; - amrex::ParticleReal p_uz_old_unit=p_uz_old[ip]; - - amrex::ParticleReal p_ux_unit=p_ux[ip]; - amrex::ParticleReal p_uy_unit=p_uy[ip]; - amrex::ParticleReal p_uz_unit=p_uz[ip]; + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip){ + amrex::ParticleReal xp, yp, zp; + GetPosition.AsStored(ip,xp, yp, zp); - amrex::ParticleReal p_u = std::sqrt(std::pow(p_ux_unit,2)+std::pow(p_uy_unit,2)+std::pow(p_uz_unit,2)); - //Calculation of gamma - amrex::ParticleReal gamma = 1/(1-std::pow(p_u,2)/std::pow(ablastr::constant::SI::c,2)); + amrex::ParticleReal u2 = p_ux[ip]*p_ux[ip]+)p_uy[ip]*p_uy[ip]+p_uz[ip]*p_uz[ip]; + amrex::ParticleReal gamma = sdt::sqrt(1.0_rt + u2); //Calculation of ei(omegat-n.r) n.resize(3); @@ -218,10 +230,17 @@ void RadiationHandler::add_radiation_contribution } }); - } + + + } } + } + + } + + 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); From 5f119731e4b505adc336ac5246be50637eb1b28f Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 19 Dec 2023 16:51:18 +0100 Subject: [PATCH 28/60] refactoring --- Source/Particles/Radiation/RadiationHandler.H | 1 - .../Particles/Radiation/RadiationHandler.cpp | 164 ++++++++++-------- 2 files changed, 96 insertions(+), 69 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 316f08c4f54..ce14f31af33 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -13,7 +13,6 @@ #include - #include #include #include diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 35b777e28f3..84eacdb0088 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -152,92 +152,120 @@ void RadiationHandler::add_radiation_contribution for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) { long const np = pti.numParticles(); - auto& attribs = pti.GetAttribs(); - auto& p_w = attribs[PIdx::w]; - auto& p_ux = attribs[PIdx::ux]; - auto& p_uy = attribs[PIdx::uy]; - auto& p_uz = attribs[PIdx::uz]; + const auto& attribs = pti.GetAttribs(); + const auto& p_w = attribs[PIdx::w]; + const auto& p_ux = attribs[PIdx::ux]; + const auto& p_uy = attribs[PIdx::uy]; + const auto& p_uz = attribs[PIdx::uz]; - auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); + const auto index = std::make_pair(pti.index(), pti.LocalTileIndex()); auto& part=pc->GetParticles(lev)[index]; auto& soa = part.GetStructOfArrays(); - auto* p_ux_old = soa.GetRealData(pc->GetRealCompIndex("old_u_x")).data(); - auto* p_uy_old = soa.GetRealData(pc->GetRealCompIndex("old_u_y")).data(); - auto* p_uz_old = soa.GetRealData(pc->GetRealCompIndex("old_u_z")).data(); + 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); - amrex::ParticleReal const q = pc->getCharge(); + 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(); + + const auto p_radiation_data = m_radiation_data.dataPtr(); 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; - amrex::ParticleReal u2 = p_ux[ip]*p_ux[ip]+)p_uy[ip]*p_uy[ip]+p_uz[ip]*p_uz[ip]; - amrex::ParticleReal gamma = sdt::sqrt(1.0_rt + u2); + constexpr auto inv_c2 = 1._prt/(PhysConst::c*PhysConst::c); + auto const one_over_gamma = 1._prt/sdt::sqrt(1.0_rt + u2*inv_c2); + auto const one_over_gamma_c = one_over_gamma/PhysConst::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; - //Calculation of ei(omegat-n.r) - n.resize(3); - //omega effective calculé - for(int i_om=0; i_om Date: Wed, 20 Dec 2023 17:05:13 +0100 Subject: [PATCH 29/60] fix issue with headers --- Source/Particles/Radiation/RadiationHandler.H | 4 ++-- Source/Particles/Radiation/RadiationHandler.cpp | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index ce14f31af33..5818ba0fb9f 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -11,6 +11,8 @@ #include "RadiationHandler_fwd.H" +#include "Particles/WarpXParticleContainer_fwd.H" + #include #include @@ -34,8 +36,6 @@ public: void gather_and_write_radiation(const std::string& filename); - void Integral_overtime(); - private: // Frequency range diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 84eacdb0088..a9728ecfce0 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -8,6 +8,7 @@ #include "RadiationHandler.H" #include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/WarpXParticleContainer.H" #include "Utils/TextMsg.H" #include From cd4db3993f42cf851336386e0388c8f60b672262 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 20 Dec 2023 17:13:14 +0100 Subject: [PATCH 30/60] added missing include --- Source/Particles/Radiation/RadiationHandler.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index a9728ecfce0..41d005dbb76 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -18,6 +18,7 @@ #include +using namespace amrex; using namespace ablastr::math; namespace From 9f6087a035caa3dba658cc71ee2e44a359e166cd Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 20 Dec 2023 17:15:57 +0100 Subject: [PATCH 31/60] fix bugs --- Source/Particles/Radiation/RadiationHandler.cpp | 4 ++-- Source/ablastr/math/Utils.H | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 41d005dbb76..f13eb8fe0cf 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -57,11 +57,11 @@ namespace auto us = amrex::Vector(det_points[0]); const auto ulim = distance*std::tan(theta_range[0]*0.5_rt); - ablastr::math::LinSpaceFill(-ul,ul, det_points[0], us.being()); + ablastr::math::LinSpaceFill(-ulim,ulim, det_points[0], us.begin()); auto vs = amrex::Vector(det_points[1]); const auto vlim = distance*std::tan(theta_range[1]*0.5_rt); - ablastr::math::LinSpaceFill(-vlim, vlim, det_points[1], vs.being()); + ablastr::math::LinSpaceFill(-vlim, vlim, det_points[1], vs.begin()); const auto how_many = det_points[0]*det_points[1]; diff --git a/Source/ablastr/math/Utils.H b/Source/ablastr/math/Utils.H index a94e1da484d..86b44786106 100644 --- a/Source/ablastr/math/Utils.H +++ b/Source/ablastr/math/Utils.H @@ -14,7 +14,7 @@ namespace ablastr::math { template AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE - void LinSpaceFill(T start, T end, int size, Iter& it) + void LinSpaceFill(T start, T end, int size, Iter it) { const auto delta = static_cast((end - start)/(size-1)); From 35003344bf4f8111a933d412a888e85dcf8a9c62 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 20 Dec 2023 17:24:08 +0100 Subject: [PATCH 32/60] fixed bugs --- .../Particles/Radiation/RadiationHandler.cpp | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index f13eb8fe0cf..6f4f23410cd 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -16,6 +16,7 @@ #include +#include #include using namespace amrex; @@ -65,9 +66,9 @@ namespace const auto how_many = det_points[0]*det_points[1]; - auto det_x = amrex::GPU::DeviceVector(how_many); - auto det_y = amrex::GPU::DeviceVector(how_many); - auto det_z = amrex::GPU::DeviceVector(how_many); + auto det_x = amrex::Gpu::DeviceVector(how_many); + auto det_y = amrex::Gpu::DeviceVector(how_many); + auto det_z = amrex::Gpu::DeviceVector(how_many); const auto one_over_direction = 1.0_rt/std::sqrt( direction[0]*direction[0]+direction[1]*direction[1]+direction[2]*direction[2]); @@ -86,9 +87,9 @@ namespace } } - amrex::GPU::synchronize(); + amrex::Gpu::synchronize(); - return {det_x, det_y, det_z} + return std::make_tuple(det_x, det_y, det_z); } } @@ -104,26 +105,26 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) const amrex::ParmParse pp_radiation("radiation"); //Resolution in frequency of the detector - auto omega_range std::vector(2); + auto omega_range = std::vector(2); pp_radiation.getarr("omega_range", omega_range); std::copy(omega_range.begin(), omega_range.end(), m_omega_range.begin()); pp_radiation.get("omega_points", m_omega_points); //Angle theta AND phi - auto theta_range std::vector(2); + auto theta_range = std::vector(2); pp_radiation.get("angle_aperture", theta_range); std::copy(theta_range.begin(), theta_range.end(), m_theta_range.begin()); //Detector parameters - auto det_pts std::vector(2); + auto det_pts = std::vector(2); pp_radiation.getarr("detector_number_points", det_pts); std::copy(det_pts.begin(), det_pts.end(), m_det_pts.begin()); - auto det_direction std::vector(3); + auto det_direction = std::vector(3); pp_radiation.getarr("detector_direction", det_direction); std::copy(det_direction.begin(), det_direction.end(), m_det_direction.begin()); - auto det_orientation std::vector(3); + auto det_orientation = std::vector(3); pp_radiation.getarr("detector_orientation", det_orientation); std::copy(det_orientation.begin(), det_orientation.end(), m_det_orientation.begin()); @@ -213,7 +214,7 @@ void RadiationHandler::add_radiation_contribution 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( + 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; From 8465c79308e8aacf1b1298227b1b0ee9275137b7 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 20 Dec 2023 17:24:39 +0100 Subject: [PATCH 33/60] fixed bugs --- Source/Particles/Radiation/RadiationHandler.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 6f4f23410cd..2c686a5399d 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -112,7 +112,7 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) //Angle theta AND phi auto theta_range = std::vector(2); - pp_radiation.get("angle_aperture", theta_range); + pp_radiation.getarr("angle_aperture", theta_range); std::copy(theta_range.begin(), theta_range.end(), m_theta_range.begin()); //Detector parameters From 2e521227bf011c23721d393f541c508664b2b4bb Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 20 Dec 2023 17:34:47 +0100 Subject: [PATCH 34/60] fixed bugs --- Source/Particles/Radiation/RadiationHandler.cpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 2c686a5399d..b3ff0b5fd58 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -16,7 +16,11 @@ #include -#include +#ifdef AMREX_USE_OMP +# include +#endif + +#include #include using namespace amrex; @@ -130,15 +134,15 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) pp_radiation.get("detector_distance", m_det_distance); - [det_pos_x, det_pos_y, det_pos_z] = compute_detector_positions( - center, det_direction, m_det_distance, + std::tie(det_pos_x, det_pos_y, det_pos_z) = compute_detector_positions( + center, m_det_direction, m_det_distance, m_det_orientation, m_det_pts, m_theta_range); constexpr auto ncomp = 3; m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); m_omegas = amrex::Gpu::DeviceVector(m_omega_points); - ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, omega_calc.being()); + ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, m_omegas.begin()); } @@ -192,7 +196,7 @@ void RadiationHandler::add_radiation_contribution auto const u2 = ux*ux + uy*uy + uz*uz; constexpr auto inv_c2 = 1._prt/(PhysConst::c*PhysConst::c); - auto const one_over_gamma = 1._prt/sdt::sqrt(1.0_rt + u2*inv_c2); + auto const one_over_gamma = 1._prt/std::sqrt(1.0_rt + u2*inv_c2); auto const one_over_gamma_c = one_over_gamma/PhysConst::c; const auto bx = ux*one_over_gamma_c; const auto by = uy*one_over_gamma_c; @@ -276,7 +280,7 @@ 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()); + m_radiation_data.begin(), m_radiation_data.end(), m_radiation_data.begin()); amrex::Gpu::streamSynchronize(); amrex::ParallelDescriptor::ReduceRealSum(radiation_data_cpu.data(), radiation_data_cpu.size()); From 58f75bcb0bf0ae13813be59330c33c21ce5f62c1 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 20 Dec 2023 17:49:30 +0100 Subject: [PATCH 35/60] fix bugs --- Source/Particles/MultiParticleContainer.cpp | 1 - .../Particles/Radiation/RadiationHandler.cpp | 25 +++++++++++++------ 2 files changed, 17 insertions(+), 9 deletions(-) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 748446525be..09d44f39f18 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -991,7 +991,6 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt, amrex::Real cur_ void MultiParticleContainer::Dump_radiations(){ if (m_at_least_one_has_radiation){ - m_p_radiation_handler->Integral_overtime(); m_p_radiation_handler->gather_and_write_radiation("Radiation"); } } diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index b3ff0b5fd58..006ebee1372 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -184,7 +184,7 @@ void RadiationHandler::add_radiation_contribution const auto p_det_pos_y = det_pos_y.dataPtr(); const auto p_det_pos_z = det_pos_z.dataPtr(); - const auto p_radiation_data = m_radiation_data.dataPtr(); + auto p_radiation_data = m_radiation_data.dataPtr(); amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip){ amrex::ParticleReal xp, yp, zp; @@ -252,18 +252,27 @@ void RadiationHandler::add_radiation_contribution const auto cz = coeff*n_cross_n_minus_beta_cross_bp_z; const int ncomp = 3; - const int idx = (i_om*how_many_det_pos + i_det)*ncomp; + 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[idx + 0] += cx; + 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[idx + 1] += cy; + p_radiation_data[idx2].m_real += cz.m_real; #pragma omp atomic - p_radiation_data[idx + 2] += cz; + p_radiation_data[idx2].m_imag += cz.m_imag; #else - p_radiation_data[idx + 0] += cx; - p_radiation_data[idx + 1] += cy; - p_radiation_data[idx + 2] += cz; + p_radiation_data[idx0] += cx; + p_radiation_data[idx1] += cy; + p_radiation_data[idx2] += cz; #endif } From e9454cdd0a6291eaf5ee3fc7e0c6dca085c4f8a2 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 11:59:45 +0100 Subject: [PATCH 36/60] fix bugs --- Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp | 3 +++ 1 file changed, 3 insertions(+) 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, From e9664b1807172424ef0f2c8a7c66ba95d22e0ad7 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 12:06:36 +0100 Subject: [PATCH 37/60] add missing line in Makefile --- Source/Particles/Radiation/Make.package | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Source/Particles/Radiation/Make.package b/Source/Particles/Radiation/Make.package index 70659a0e305..1c3a1020746 100644 --- a/Source/Particles/Radiation/Make.package +++ b/Source/Particles/Radiation/Make.package @@ -1 +1,3 @@ CEXE_sources += RadiationHandler.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Radiation From 3d6cc2161101f0dabd5f4ee72ea7c5d6b3da2b66 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 12:14:41 +0100 Subject: [PATCH 38/60] fix misplaced using namespace directive --- Source/Particles/Deposition/SharedDepositionUtils.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Deposition/SharedDepositionUtils.H b/Source/Particles/Deposition/SharedDepositionUtils.H index f3da11e94a9..cb6a9e4cdbe 100644 --- a/Source/Particles/Deposition/SharedDepositionUtils.H +++ b/Source/Particles/Deposition/SharedDepositionUtils.H @@ -84,9 +84,9 @@ void depositComponent (const GetParticlePosition& GetPosition, const int zdir, const int NODE, const int CELL, const int dir) { using namespace amrex::literals; - using namespace ablastr::math; #if !defined(WARPX_DIM_RZ) + using namespace ablastr::math; amrex::ignore_unused(n_rz_azimuthal_modes); #endif From 8572d2b2147fcd38addda59e2e04a57dc0d972dc Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 13:57:33 +0100 Subject: [PATCH 39/60] add missing include --- Source/FieldSolver/SpectralSolver/SpectralKSpace.cpp | 3 +++ 1 file changed, 3 insertions(+) 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. * From 48fc5a7274178736e8dbb2ff761850b8c92bc203 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 14:07:24 +0100 Subject: [PATCH 40/60] fix misspelled name --- .../SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp index 3e564f74b31..c00ab2ff79e 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp @@ -28,7 +28,7 @@ #if WARPX_USE_PSATD using namespace amrex; -using namespace ablsatr::math; +using namespace ablastr::math; PsatdAlgorithmJConstantInTime::PsatdAlgorithmJConstantInTime( const SpectralKSpace& spectral_kspace, From e42a20aa50dcd9cc4afe9b066df0545c4950b22f Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 14:16:31 +0100 Subject: [PATCH 41/60] fixed bug --- Source/Particles/Deposition/SharedDepositionUtils.H | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/Source/Particles/Deposition/SharedDepositionUtils.H b/Source/Particles/Deposition/SharedDepositionUtils.H index cb6a9e4cdbe..3b46c340f35 100644 --- a/Source/Particles/Deposition/SharedDepositionUtils.H +++ b/Source/Particles/Deposition/SharedDepositionUtils.H @@ -86,8 +86,9 @@ void depositComponent (const GetParticlePosition& GetPosition, using namespace amrex::literals; #if !defined(WARPX_DIM_RZ) - using namespace ablastr::math; 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 From 7475fd8a53721182b937a70a218745335c421a07 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 14:45:37 +0100 Subject: [PATCH 42/60] using parser to read inputfile --- .../Particles/Radiation/RadiationHandler.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 006ebee1372..6f5e1247718 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -9,13 +9,12 @@ #include "Particles/Pusher/GetAndSetPosition.H" #include "Particles/WarpXParticleContainer.H" +#include "Utils/Parser/ParserUtils.H" #include "Utils/TextMsg.H" #include #include -#include - #ifdef AMREX_USE_OMP # include #endif @@ -25,6 +24,7 @@ using namespace amrex; using namespace ablastr::math; +using namespace utils::parser; namespace { @@ -110,29 +110,29 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) //Resolution in frequency of the detector auto omega_range = std::vector(2); - pp_radiation.getarr("omega_range", omega_range); + getArrWithParser(pp_radiation, "omega_range", omega_range); std::copy(omega_range.begin(), omega_range.end(), m_omega_range.begin()); - pp_radiation.get("omega_points", m_omega_points); + getWithParser(pp_radiation, "omega_points", m_omega_points); //Angle theta AND phi auto theta_range = std::vector(2); - pp_radiation.getarr("angle_aperture", theta_range); + 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); - pp_radiation.getarr("detector_number_points", det_pts); + 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); - pp_radiation.getarr("detector_direction", det_direction); + 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); - pp_radiation.getarr("detector_orientation", det_orientation); + getArrWithParser(pp_radiation, "detector_orientation", det_orientation); std::copy(det_orientation.begin(), det_orientation.end(), m_det_orientation.begin()); - pp_radiation.get("detector_distance", m_det_distance); + getWithParser(pp_radiation, "detector_distance", m_det_distance); std::tie(det_pos_x, det_pos_y, det_pos_z) = compute_detector_positions( center, m_det_direction, m_det_distance, From cbcef567f92673a38fef0e2855d02ac8a754cabe Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 14:47:14 +0100 Subject: [PATCH 43/60] fix bug --- Source/Particles/Radiation/RadiationHandler.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 6f5e1247718..a7478a6dec7 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -86,8 +86,8 @@ namespace for (int j = 0; j < det_points[1]; ++j) { det_x[i*det_points[1] + j] = center[0] + distance * norm_direction[0] + us[i]*u[0] + vs[j]*v[0]; - det_x[i*det_points[1] + j] = center[1] + distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; - det_x[i*det_points[1] + j] = center[2] + distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; + det_y[i*det_points[1] + j] = center[1] + distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; + det_z[i*det_points[1] + j] = center[2] + distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; } } From 7b94a13e22d910eb84d480ccdd6a35e40ff5032b Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:18:22 +0100 Subject: [PATCH 44/60] make function public to allow its use on NVIDIA GPUs --- Source/Particles/WarpXParticleContainer.H | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 7f166c3d3d1..71c202a342f 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -406,6 +406,12 @@ 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: /** @@ -414,11 +420,6 @@ protected: */ void AllocateOldMomenta(); - /** - * This function copies the momenta of the particles into: - * old_u_x, old_u_y, old_u_z. - */ - void RecordOldMomenta(); int species_id; From 53ea810619bc92bdfb4a03a582bb57c41c8e9794 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:23:40 +0100 Subject: [PATCH 45/60] avoid implicit capture of this --- Source/Particles/Radiation/RadiationHandler.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index a7478a6dec7..5c732aafe32 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -186,6 +186,8 @@ void RadiationHandler::add_radiation_contribution auto p_radiation_data = m_radiation_data.dataPtr(); + const auto& omega_points = m_omega_points; + amrex::ParallelFor(np, [=] AMREX_GPU_DEVICE(int ip){ amrex::ParticleReal xp, yp, zp; GetPosition.AsStored(ip,xp, yp, zp); @@ -209,7 +211,7 @@ void RadiationHandler::add_radiation_contribution const auto tot_q = q*p_w[ip]; - for(int i_om=0; i_om < m_omega_points; ++i_om){ + 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]/PhysConst::c; From 6af40a09da5a045f3a754227994e2e4c3b08f926 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:26:45 +0100 Subject: [PATCH 46/60] make LinSpaceFill callable only from host code --- Source/ablastr/math/Utils.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/ablastr/math/Utils.H b/Source/ablastr/math/Utils.H index 86b44786106..88f4c195a62 100644 --- a/Source/ablastr/math/Utils.H +++ b/Source/ablastr/math/Utils.H @@ -13,7 +13,7 @@ namespace ablastr::math { template - AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE + AMREX_FORCE_INLINE void LinSpaceFill(T start, T end, int size, Iter it) { const auto delta = static_cast((end - start)/(size-1)); From ab65f9efa3d6971f810ed6ceeb5e31d9113ecfa5 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:29:20 +0100 Subject: [PATCH 47/60] use pointers for data to be used on GPUs --- Source/Particles/Radiation/RadiationHandler.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 5c732aafe32..bc31dfa6133 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -160,10 +160,10 @@ void RadiationHandler::add_radiation_contribution { long const np = pti.numParticles(); const auto& attribs = pti.GetAttribs(); - const auto& p_w = attribs[PIdx::w]; - const auto& p_ux = attribs[PIdx::ux]; - const auto& p_uy = attribs[PIdx::uy]; - const auto& p_uz = attribs[PIdx::uz]; + 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]; From 753409fb996d970b5abe93353335f0c704e9c487 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:33:20 +0100 Subject: [PATCH 48/60] missing include --- Source/Particles/Radiation/RadiationHandler.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index bc31dfa6133..ffd1dda634e 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -10,6 +10,7 @@ #include "Particles/Pusher/GetAndSetPosition.H" #include "Particles/WarpXParticleContainer.H" #include "Utils/Parser/ParserUtils.H" +#include "Utils/WarpXConst.H" #include "Utils/TextMsg.H" #include From 4faf90ed0c585f675b98c2685590514f095b265e Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:46:40 +0100 Subject: [PATCH 49/60] removing include --- Source/Particles/Radiation/RadiationHandler.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index ffd1dda634e..155a58d908c 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -13,7 +13,6 @@ #include "Utils/WarpXConst.H" #include "Utils/TextMsg.H" -#include #include #ifdef AMREX_USE_OMP From 9c25bbf0048c20f6423912ef299a70c44f945dce Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 15:54:50 +0100 Subject: [PATCH 50/60] fix issue with constexpr --- Source/Particles/Radiation/RadiationHandler.cpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 155a58d908c..50ea5cf3fe5 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -188,6 +188,10 @@ void RadiationHandler::add_radiation_contribution 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); @@ -197,9 +201,8 @@ void RadiationHandler::add_radiation_contribution const auto uz = 0.5_prt*(p_uz[ip] + p_uz_old[ip]); auto const u2 = ux*ux + uy*uy + uz*uz; - constexpr auto inv_c2 = 1._prt/(PhysConst::c*PhysConst::c); auto const one_over_gamma = 1._prt/std::sqrt(1.0_rt + u2*inv_c2); - auto const one_over_gamma_c = one_over_gamma/PhysConst::c; + 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; @@ -213,7 +216,7 @@ void RadiationHandler::add_radiation_contribution 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]/PhysConst::c; + 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){ @@ -245,7 +248,7 @@ void RadiationHandler::add_radiation_contribution 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*(PhysConst::c*current_time - d_part_det)); + 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); From e53cb0891de8c6eb81dc23fd52f9e5a2a57318ec Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 21 Dec 2023 16:11:01 +0100 Subject: [PATCH 51/60] fix issues on GPU --- .../Particles/Radiation/RadiationHandler.cpp | 34 ++++++++++++++----- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 50ea5cf3fe5..bb9c897f7e9 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -70,9 +70,9 @@ namespace const auto how_many = det_points[0]*det_points[1]; - auto det_x = amrex::Gpu::DeviceVector(how_many); - auto det_y = amrex::Gpu::DeviceVector(how_many); - auto det_z = amrex::Gpu::DeviceVector(how_many); + auto host_det_x = amrex::Vector(how_many); + auto host_det_y = amrex::Vector(how_many); + auto host_det_z = amrex::Vector(how_many); const auto one_over_direction = 1.0_rt/std::sqrt( direction[0]*direction[0]+direction[1]*direction[1]+direction[2]*direction[2]); @@ -85,15 +85,26 @@ namespace { for (int j = 0; j < det_points[1]; ++j) { - det_x[i*det_points[1] + j] = center[0] + distance * norm_direction[0] + us[i]*u[0] + vs[j]*v[0]; - det_y[i*det_points[1] + j] = center[1] + distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; - det_z[i*det_points[1] + j] = center[2] + distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; + host_det_x[i*det_points[1] + j] = center[0] + distance * norm_direction[0] + us[i]*u[0] + vs[j]*v[0]; + host_det_y[i*det_points[1] + j] = center[1] + distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; + host_det_z[i*det_points[1] + j] = center[2] + distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; } } - amrex::Gpu::synchronize(); + 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); - return std::make_tuple(det_x, det_y, det_z); + 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::Device::streamSynchronize(); + + return std::make_tuple(gpu_det_x, gpu_det_y, gpu_det_z); } } @@ -141,8 +152,13 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) 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); + ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, t_omegas.begin()); m_omegas = amrex::Gpu::DeviceVector(m_omega_points); - ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, m_omegas.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, + t_omegas.begin(), t_omegas.end(), m_omegas.begin()); + + amrex::Gpu::Device::streamSynchronize(); } From 06a5a111fe396b6b7f46677f43ff1345fdd33418 Mon Sep 17 00:00:00 2001 From: tmsclark Date: Tue, 23 Jan 2024 14:08:13 +0100 Subject: [PATCH 52/60] Integral_over_time --- Source/Evolve/WarpXEvolve.cpp | 3 +- Source/Particles/MultiParticleContainer.H | 2 +- Source/Particles/MultiParticleContainer.cpp | 4 +- Source/Particles/Radiation/RadiationHandler.H | 4 ++ .../Particles/Radiation/RadiationHandler.cpp | 37 ++++++++++++------- 5 files changed, 33 insertions(+), 17 deletions(-) diff --git a/Source/Evolve/WarpXEvolve.cpp b/Source/Evolve/WarpXEvolve.cpp index f7cef17c8e2..e57abea36ce 100644 --- a/Source/Evolve/WarpXEvolve.cpp +++ b/Source/Evolve/WarpXEvolve.cpp @@ -246,7 +246,8 @@ WarpX::Evolve (int numsteps) *Bfield_aux[lev][0],*Bfield_aux[lev][1], *Bfield_aux[lev][2]); } - mypc->Dump_radiations(); + // For now only on lower level + mypc->Dump_radiations(dt[0]); is_synchronized = true; } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 85b95894513..791f7c77329 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -334,7 +334,7 @@ public: void RecordOldMomenta(); - void Dump_radiations(); + void Dump_radiations(amrex::Real dt); [[nodiscard]] int getSpeciesID (std::string product_str) const; protected: diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 09d44f39f18..7970511dbfd 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -988,9 +988,9 @@ void MultiParticleContainer::doRadiation (const amrex::Real dt, amrex::Real cur_ } } - -void MultiParticleContainer::Dump_radiations(){ +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"); } } diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 5818ba0fb9f..ea67bbb0c47 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -36,6 +36,8 @@ public: void gather_and_write_radiation(const std::string& filename); + void Integral_overtime(const amrex::Real dt); + private: // Frequency range @@ -56,6 +58,8 @@ private: amrex::Gpu::DeviceVector det_pos_z; 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 index bb9c897f7e9..e00306f8b39 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -79,7 +79,7 @@ namespace const auto norm_direction = amrex::Array{ direction[0]*one_over_direction, direction[1]*one_over_direction, - direction[2]*one_over_direction}; + direction[2]*one_over_direction}; for (int i = 0; i < det_points[0]; ++i) { @@ -159,19 +159,21 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) 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 - { + { amrex::Print() << "Je brille <<" << std::endl; + for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) { long const np = pti.numParticles(); @@ -231,7 +233,7 @@ void RadiationHandler::add_radiation_contribution 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){ @@ -258,6 +260,7 @@ void RadiationHandler::add_radiation_contribution 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; + amrex::Print() << "A" << bx << "b" << p_ux_old << std::endl; //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; @@ -294,7 +297,7 @@ void RadiationHandler::add_radiation_contribution p_radiation_data[idx0] += cx; p_radiation_data[idx1] += cy; p_radiation_data[idx2] += cz; - + #endif } } @@ -310,7 +313,7 @@ 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_data.begin(), m_radiation_data.end(), m_radiation_data.begin()); + m_radiation_calculation.begin(), m_radiation_calculation.end(), m_radiation_calculation.begin()); amrex::Gpu::streamSynchronize(); amrex::ParallelDescriptor::ReduceRealSum(radiation_data_cpu.data(), radiation_data_cpu.size()); @@ -326,13 +329,21 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) } } -/*void RadiationHandler::Integral_overtime() -{ +void RadiationHandler::Integral_overtime(const amrex::Real dt) +{ + const auto factor = dt/16/std::pow(ablastr::constant::math::pi,3)/PhysConst::ep0; + const auto how_many = m_det_pts[0]*m_det_pts[1]; + auto p_radiation_data = m_radiation_data.dataPtr(); for(int i_om=0; i_om Date: Tue, 23 Jan 2024 14:25:36 +0100 Subject: [PATCH 53/60] Integral_over_time_mod --- Source/Particles/Radiation/RadiationHandler.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index e00306f8b39..df17a8d96d0 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -260,7 +260,6 @@ void RadiationHandler::add_radiation_contribution 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; - amrex::Print() << "A" << bx << "b" << p_ux_old << std::endl; //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; From 039957c94d977cbd6b3d15f1fc31775d1819a153 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 23 Jan 2024 15:53:56 +0100 Subject: [PATCH 54/60] debugging radiation module --- .../Particles/Radiation/RadiationHandler.cpp | 48 ++++++++++++------- 1 file changed, 32 insertions(+), 16 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index df17a8d96d0..ea4200c16ef 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -39,7 +39,7 @@ namespace WARPX_ALWAYS_ASSERT_WITH_MESSAGE( direction[0]*orientation[0] + direction[1]*orientation[1] + - direction[2]*orientation[2] != 0, + direction[2]*orientation[2] == 0, "Radiation detector orientation cannot be aligned with detector direction"); auto u = amrex::Array{ @@ -79,7 +79,7 @@ namespace const auto norm_direction = amrex::Array{ direction[0]*one_over_direction, direction[1]*one_over_direction, - direction[2]*one_over_direction}; + direction[2]*one_over_direction}; for (int i = 0; i < det_points[0]; ++i) { @@ -233,7 +233,7 @@ void RadiationHandler::add_radiation_contribution 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){ @@ -296,7 +296,6 @@ void RadiationHandler::add_radiation_contribution p_radiation_data[idx0] += cx; p_radiation_data[idx1] += cy; p_radiation_data[idx2] += cz; - #endif } } @@ -312,7 +311,7 @@ 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(), m_radiation_calculation.begin()); + 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()); @@ -320,8 +319,29 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) if (amrex::ParallelDescriptor::IOProcessor() ){ auto of = std::ofstream(filename, std::ios::binary); - for (const auto& d : radiation_data_cpu){ - of << static_cast(d); + 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 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, + 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_x_cpu[i_det] << " " << det_pos_y_cpu[i_det] << " " << det_pos_z_cpu[i_det] << " " << radiation_data_cpu[++idx] << "\n"; + } } of.close(); @@ -329,20 +349,16 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) } void RadiationHandler::Integral_overtime(const amrex::Real dt) -{ +{ const auto factor = dt/16/std::pow(ablastr::constant::math::pi,3)/PhysConst::ep0; const auto how_many = m_det_pts[0]*m_det_pts[1]; auto p_radiation_data = m_radiation_data.dataPtr(); - for(int i_om=0; i_om Date: Wed, 31 Jan 2024 11:00:05 +0100 Subject: [PATCH 55/60] Nyquist limiter --- .../Particles/Radiation/RadiationHandler.cpp | 22 +++++++++++-------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index ea4200c16ef..f4be3aca9f0 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -172,8 +172,7 @@ void RadiationHandler::add_radiation_contribution #ifdef AMREX_USE_OMP #pragma omp parallel #endif - { amrex::Print() << "Je brille <<" << std::endl; - + { for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) { long const np = pti.numParticles(); @@ -270,15 +269,21 @@ void RadiationHandler::add_radiation_contribution const auto coeff = tot_q*phase_term/(one_minus_b_dot_n*one_minus_b_dot_n); - const auto cx = coeff*n_cross_n_minus_beta_cross_bp_x; - const auto cy = coeff*n_cross_n_minus_beta_cross_bp_y; - const auto cz = coeff*n_cross_n_minus_beta_cross_bp_z; - 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; + 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(ablastr::constant::math::pi/(dt*one_minus_b_dot_n)*p_omegas[i_om] < 0){ + cx = 0.0; + cy = 0.0; + cz = 0.0; + } #ifdef AMREX_USE_OMP #pragma omp atomic p_radiation_data[idx0].m_real += cx.m_real; @@ -294,7 +299,7 @@ void RadiationHandler::add_radiation_contribution p_radiation_data[idx2].m_imag += cz.m_imag; #else p_radiation_data[idx0] += cx; - p_radiation_data[idx1] += cy; + p_radiation_data[idx1 ] += cy; p_radiation_data[idx2] += cz; #endif } @@ -350,7 +355,7 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) void RadiationHandler::Integral_overtime(const amrex::Real dt) { - const auto factor = dt/16/std::pow(ablastr::constant::math::pi,3)/PhysConst::ep0; + const auto factor = 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); @@ -358,7 +363,6 @@ void RadiationHandler::Integral_overtime(const amrex::Real dt) const int idx0 = idx*3; const int idx1 = idx0 + 1; const int idx2 = idx0 + 2; - amrex::Print() << (amrex::norm(p_radiation_data[idx0])+amrex::norm(p_radiation_data[idx1])+amrex::norm(p_radiation_data[idx2])) << std::endl; m_radiation_calculation[idx]=(amrex::norm(p_radiation_data[idx0])+amrex::norm(p_radiation_data[idx1])+amrex::norm(p_radiation_data[idx2]))*factor; } } From f3d4a295a65124e90d2914d9a520cf7c4b559ea9 Mon Sep 17 00:00:00 2001 From: Pierre Bartoli Date: Wed, 7 Feb 2024 10:25:46 +0100 Subject: [PATCH 56/60] Work in progress: spherical detector --- Source/Particles/Radiation/RadiationHandler.H | 2 + .../Particles/Radiation/RadiationHandler.cpp | 119 ++++++++++++++---- 2 files changed, 99 insertions(+), 22 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index ea67bbb0c47..205ac829d20 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -56,6 +56,8 @@ private: 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; diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index f4be3aca9f0..3866dccf77d 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -36,7 +36,7 @@ namespace const amrex::Array& det_points, const amrex::Array& theta_range) { - WARPX_ALWAYS_ASSERT_WITH_MESSAGE( + /*WARPX_ALWAYS_ASSERT_WITH_MESSAGE( direction[0]*orientation[0] + direction[1]*orientation[1] + direction[2]*orientation[2] == 0, @@ -105,6 +105,71 @@ namespace amrex::Gpu::Device::streamSynchronize(); return std::make_tuple(gpu_det_x, gpu_det_y, gpu_det_z); + */ + + 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[0]) - Ntheta/2*dtheta; + const auto phiOrigin = std::atan2(direction[1], direction[0]) - Nphi/2*dphi; + + 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); + + 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::sin(phi); + host_det_y[i*det_points[1] + j] = center[1] + distance*std::sin(theta)*std::cos(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; + + //amrex::Print() << phi << std::endl; + } + } + + 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); } } @@ -145,7 +210,7 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) getWithParser(pp_radiation, "detector_distance", m_det_distance); - std::tie(det_pos_x, det_pos_y, det_pos_z) = compute_detector_positions( + 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); @@ -157,7 +222,6 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) 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(); } @@ -172,7 +236,8 @@ void RadiationHandler::add_radiation_contribution #ifdef AMREX_USE_OMP #pragma omp parallel #endif - { + { + for (WarpXParIter pti(*pc, lev); pti.isValid(); ++pti) { long const np = pti.numParticles(); @@ -207,7 +272,7 @@ void RadiationHandler::add_radiation_contribution constexpr auto c = PhysConst::c; constexpr auto inv_c = 1._prt/(PhysConst::c); - constexpr auto inv_c2 = 1._prt/(PhysConst::c* 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; @@ -219,6 +284,7 @@ void RadiationHandler::add_radiation_contribution 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; @@ -231,10 +297,11 @@ void RadiationHandler::add_radiation_contribution const auto tot_q = q*p_w[ip]; + //amrex::Print() << tot_q << std::endl; + 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]; @@ -261,19 +328,14 @@ void RadiationHandler::add_radiation_contribution 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 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); - 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; - 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; @@ -283,7 +345,13 @@ void RadiationHandler::add_radiation_contribution cx = 0.0; cy = 0.0; cz = 0.0; + amrex::Print() << cx << std::endl; } + 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; @@ -299,7 +367,7 @@ void RadiationHandler::add_radiation_contribution p_radiation_data[idx2].m_imag += cz.m_imag; #else p_radiation_data[idx0] += cx; - p_radiation_data[idx1 ] += cy; + p_radiation_data[idx1] += cy; p_radiation_data[idx2] += cz; #endif } @@ -325,17 +393,23 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) 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_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_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, + /*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()); + 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(); @@ -345,7 +419,7 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) 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_x_cpu[i_det] << " " << det_pos_y_cpu[i_det] << " " << det_pos_z_cpu[i_det] << " " << radiation_data_cpu[++idx] << "\n"; + of << omegas_cpu[i_om] << " " << det_pos_theta_cpu[i_det] << " " << det_pos_phi_cpu[i_det] << " " << radiation_data_cpu[++idx] << "\n"; } } @@ -355,7 +429,7 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) void RadiationHandler::Integral_overtime(const amrex::Real dt) { - const auto factor = dt/16/std::pow(ablastr::constant::math::pi,3)/PhysConst::ep0/PhysConst::c; + const auto factor = 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); @@ -363,6 +437,7 @@ void RadiationHandler::Integral_overtime(const amrex::Real dt) const int idx0 = idx*3; const int idx1 = idx0 + 1; const int idx2 = idx0 + 2; - m_radiation_calculation[idx]=(amrex::norm(p_radiation_data[idx0])+amrex::norm(p_radiation_data[idx1])+amrex::norm(p_radiation_data[idx2]))*factor; + //amrex::Print() << (amrex::norm(p_radiation_data[idx0])+amrex::norm(p_radiation_data[idx1])+amrex::norm(p_radiation_data[idx2])) << std::endl; + m_radiation_calculation[idx]=(amrex::norm(p_radiation_data[idx0]) + amrex::norm(p_radiation_data[idx1]) + amrex::norm(p_radiation_data[idx2]))*factor; } } From 5359092bda235b4c9ecb2a2093ca7059e06a0b0a Mon Sep 17 00:00:00 2001 From: Pierre Bartoli Date: Wed, 7 Feb 2024 13:55:27 +0100 Subject: [PATCH 57/60] Correction of spherical radiation and chose between spherical and cartesian --- Source/Particles/Radiation/RadiationHandler.H | 1 - .../Particles/Radiation/RadiationHandler.cpp | 201 +++++++++--------- 2 files changed, 97 insertions(+), 105 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.H b/Source/Particles/Radiation/RadiationHandler.H index 205ac829d20..1b1a2c8f3a0 100644 --- a/Source/Particles/Radiation/RadiationHandler.H +++ b/Source/Particles/Radiation/RadiationHandler.H @@ -62,6 +62,5 @@ private: 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 index 3866dccf77d..7f3c5489215 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -34,118 +34,107 @@ namespace const amrex::Real distance, const amrex::Array& orientation, const amrex::Array& det_points, - const amrex::Array& theta_range) + const amrex::Array& theta_range, + const std::string type) { - /*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); - ablastr::math::LinSpaceFill(-ulim,ulim, det_points[0], us.begin()); - - auto vs = amrex::Vector(det_points[1]); - const auto vlim = distance*std::tan(theta_range[1]*0.5_rt); - ablastr::math::LinSpaceFill(-vlim, vlim, det_points[1], vs.begin()); 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); - - 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) - { - host_det_x[i*det_points[1] + j] = center[0] + distance * norm_direction[0] + us[i]*u[0] + vs[j]*v[0]; - host_det_y[i*det_points[1] + j] = center[1] + distance * norm_direction[1] + us[i]*u[1] + vs[j]*v[1]; - host_det_z[i*det_points[1] + j] = center[2] + distance * norm_direction[2] + us[i]*u[2] + vs[j]*v[2]; - } - } - - 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); - - 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::Device::streamSynchronize(); - - return std::make_tuple(gpu_det_x, gpu_det_y, gpu_det_z); - */ - - 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[0]) - Ntheta/2*dtheta; - const auto phiOrigin = std::atan2(direction[1], direction[0]) - Nphi/2*dphi; - - 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); - for (int i = 0; i < det_points[0]; ++i) - { - for (int j = 0; j < det_points[1]; ++j) + if(type == "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); + ablastr::math::LinSpaceFill(-ulim,ulim, det_points[0], us.begin()); + + auto vs = amrex::Vector(det_points[1]); + const auto vlim = distance*std::tan(theta_range[1]*0.5_rt); + ablastr::math::LinSpaceFill(-vlim, vlim, det_points[1], vs.begin()); + + 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) { - 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)); + 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_x[i*det_points[1] + j] = center[0] + distance*std::sin(theta)*std::sin(phi); - host_det_y[i*det_points[1] + j] = center[1] + distance*std::sin(theta)*std::cos(phi); - host_det_z[i*det_points[1] + j] = center[2] + distance*std::cos(theta); + host_det_theta[i*det_points[1] + j] = std::acos(z/distance); + host_det_phi[i*det_points[1] + j] = std::atan2(y, x); - host_det_theta[i*det_points[1] + j] = theta; - host_det_phi[i*det_points[1] + j] = phi; - - //amrex::Print() << phi << std::endl; + } + } + } + + if(type == "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; + } } } @@ -184,6 +173,10 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) // Read in radiation input const amrex::ParmParse pp_radiation("radiation"); + //type of detector + std::string type = "spherical"; + pp_radiation.query("detector_type", type); + //Resolution in frequency of the detector auto omega_range = std::vector(2); getArrWithParser(pp_radiation, "omega_range", omega_range); @@ -212,7 +205,7 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) 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); + m_det_orientation, m_det_pts, m_theta_range, type); constexpr auto ncomp = 3; m_radiation_data = amrex::Gpu::DeviceVector(m_det_pts[0]*m_det_pts[1]*m_omega_points*ncomp); @@ -393,19 +386,19 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) 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_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_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, + 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());*/ + 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, @@ -419,7 +412,7 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) 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] << " " << radiation_data_cpu[++idx] << "\n"; + 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"; } } From 9def939b176f41747791d651fad5940780d9e1a8 Mon Sep 17 00:00:00 2001 From: Pierre Bartoli Date: Mon, 12 Feb 2024 15:57:09 +0100 Subject: [PATCH 58/60] Corrections, need change of frame in spherical implementation --- Source/Particles/Radiation/RadiationHandler.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index 7f3c5489215..b7b5d052b66 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -47,6 +47,7 @@ namespace auto host_det_phi = amrex::Vector(how_many); if(type == "cartesian"){ + WARPX_ALWAYS_ASSERT_WITH_MESSAGE( direction[0]*orientation[0] + direction[1]*orientation[1] + @@ -99,13 +100,14 @@ namespace 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); + host_det_phi[i*det_points[1] + j] = std::atan2(y, x) + ablastr::constant::math::pi; } } } if(type == "spherical"){ + const auto Ntheta = det_points[0]; const auto Nphi = det_points[1]; const auto dtheta = theta_range[0]/Ntheta; @@ -290,8 +292,6 @@ void RadiationHandler::add_radiation_contribution const auto tot_q = q*p_w[ip]; - //amrex::Print() << tot_q << std::endl; - 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; @@ -334,12 +334,12 @@ void RadiationHandler::add_radiation_contribution auto cz = coeff*n_cross_n_minus_beta_cross_bp_z; // Nyquist limiter - if(ablastr::constant::math::pi/(dt*one_minus_b_dot_n)*p_omegas[i_om] < 0){ + /*if(p_omegas[i_om] < ablastr::constant::math::pi/one_minus_b_dot_n/dt){ cx = 0.0; cy = 0.0; cz = 0.0; - amrex::Print() << cx << std::endl; - } + }*/ + const int ncomp = 3; const int idx0 = (i_om*how_many_det_pos + i_det)*ncomp; const int idx1 = idx0 + 1; @@ -422,7 +422,7 @@ void RadiationHandler::gather_and_write_radiation(const std::string& filename) void RadiationHandler::Integral_overtime(const amrex::Real dt) { - const auto factor = dt/16/std::pow(ablastr::constant::math::pi,3)/PhysConst::ep0/(PhysConst::c); + 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); @@ -432,5 +432,5 @@ void RadiationHandler::Integral_overtime(const amrex::Real dt) const int idx2 = idx0 + 2; //amrex::Print() << (amrex::norm(p_radiation_data[idx0])+amrex::norm(p_radiation_data[idx1])+amrex::norm(p_radiation_data[idx2])) << std::endl; m_radiation_calculation[idx]=(amrex::norm(p_radiation_data[idx0]) + amrex::norm(p_radiation_data[idx1]) + amrex::norm(p_radiation_data[idx2]))*factor; - } + } } From 08dd0d51915edfd59c82d57c7ea5a7019d2611d4 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 15 Feb 2024 11:32:03 +0100 Subject: [PATCH 59/60] use amrex::linspace --- .../Particles/Radiation/RadiationHandler.cpp | 17 +++++------ Source/ablastr/math/Utils.H | 30 ------------------- cmake/dependencies/AMReX.cmake | 2 +- 3 files changed, 8 insertions(+), 41 deletions(-) delete mode 100644 Source/ablastr/math/Utils.H diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index b7b5d052b66..965d26ae977 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -13,7 +13,7 @@ #include "Utils/WarpXConst.H" #include "Utils/TextMsg.H" -#include +#include #ifdef AMREX_USE_OMP # include @@ -74,11 +74,11 @@ namespace auto us = amrex::Vector(det_points[0]); const auto ulim = distance*std::tan(theta_range[0]*0.5_rt); - ablastr::math::LinSpaceFill(-ulim,ulim, det_points[0], us.begin()); + 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); - ablastr::math::LinSpaceFill(-vlim, vlim, det_points[1], vs.begin()); + 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]); @@ -105,7 +105,7 @@ namespace } } } - + if(type == "spherical"){ const auto Ntheta = det_points[0]; @@ -121,7 +121,7 @@ namespace { 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; @@ -213,25 +213,23 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) 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); - ablastr::math::LinSpaceFill(m_omega_range[0], m_omega_range[1], m_omega_points, t_omegas.begin()); + 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) { @@ -339,7 +337,6 @@ void RadiationHandler::add_radiation_contribution 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; diff --git a/Source/ablastr/math/Utils.H b/Source/ablastr/math/Utils.H deleted file mode 100644 index 88f4c195a62..00000000000 --- a/Source/ablastr/math/Utils.H +++ /dev/null @@ -1,30 +0,0 @@ -/* Copyright 2023 Thomas Clark, Luca Fedeli - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#ifndef ABLASTR_MATH_UTILS_H_ -#define ABLASTR_MATH_UTILS_H_ - -#include -#include - -namespace ablastr::math -{ - template - AMREX_FORCE_INLINE - void LinSpaceFill(T start, T end, int size, Iter it) - { - const auto delta = static_cast((end - start)/(size-1)); - - for (int i = 0; i < size-1; ++i) - { - *(it++) = start + static_cast(i*delta); - } - *it = end; - } - -} - -#endif //ABLASTR_MATH_UTILS_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)") From 43d82ecc586e0975b883fc7f2c276dd47f3fbc4f Mon Sep 17 00:00:00 2001 From: Pierre Bartoli Date: Thu, 15 Feb 2024 17:27:59 +0100 Subject: [PATCH 60/60] Minor correction and adjustements --- .../Particles/Radiation/RadiationHandler.cpp | 27 +++++++++++++------ 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/Source/Particles/Radiation/RadiationHandler.cpp b/Source/Particles/Radiation/RadiationHandler.cpp index b7b5d052b66..9aac5bc329d 100644 --- a/Source/Particles/Radiation/RadiationHandler.cpp +++ b/Source/Particles/Radiation/RadiationHandler.cpp @@ -26,6 +26,16 @@ 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( @@ -35,7 +45,7 @@ namespace const amrex::Array& orientation, const amrex::Array& det_points, const amrex::Array& theta_range, - const std::string type) + const std::string radiation_type) { const auto how_many = det_points[0]*det_points[1]; @@ -46,7 +56,7 @@ namespace auto host_det_theta = amrex::Vector(how_many); auto host_det_phi = amrex::Vector(how_many); - if(type == "cartesian"){ + if(m_radiation_type.at(radiation_type) == RadiationType::cartesian){ WARPX_ALWAYS_ASSERT_WITH_MESSAGE( direction[0]*orientation[0] + @@ -106,7 +116,7 @@ namespace } } - if(type == "spherical"){ + if(m_radiation_type.at(radiation_type) == RadiationType::spherical){ const auto Ntheta = det_points[0]; const auto Nphi = det_points[1]; @@ -176,8 +186,9 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) const amrex::ParmParse pp_radiation("radiation"); //type of detector - std::string type = "spherical"; - pp_radiation.query("detector_type", type); + 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); @@ -207,7 +218,7 @@ RadiationHandler::RadiationHandler(const amrex::Array& center) 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, type); + 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); @@ -243,7 +254,7 @@ void RadiationHandler::add_radiation_contribution 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& part = pc->GetParticles(lev)[index]; auto& soa = part.GetStructOfArrays(); const auto* p_ux_old = soa.GetRealData(pc->GetRealCompIndex("old_u_x")).data(); @@ -339,7 +350,7 @@ void RadiationHandler::add_radiation_contribution 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;