From 2282ebda111ff12b8eb00c8d9c2099c20a0db99e Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Mon, 18 Dec 2023 16:04:51 +0100 Subject: [PATCH 1/3] 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 2/3] 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 3/3] 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