[WIP] Add a Radiation module

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 <ablastr/utils/SignalHandling.H>
 #include <ablastr/warn_manager/WarnManager.H>
@@ -245,6 +246,7 @@ WarpX::Evolve (int numsteps)
                             *Bfield_aux[lev][0],*Bfield_aux[lev][1],
                             *Bfield_aux[lev][2]);
             }
+            mypc->Dump_radiations();
             is_synchronized = true;
         }
 
@@ -426,8 +428,20 @@ WarpX::OneStep_nosub (Real cur_time)
     ExecutePythonCallback("particlescraper");
     ExecutePythonCallback("beforedeposition");
 
+    //Save particle old momentum in a attribute
+    mypc->keepoldmomentum();
+
     PushParticlesandDeposit(cur_time);
 
+    //Radiation contribution at each timestep
+   //Only level 0 is supported
+    mypc->doRadiation(dt[0],cur_time);
+
+
+
+
+
+
     ExecutePythonCallback("afterdeposition");
 
     // Synchronize J and rho:

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)

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
 

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,12 @@ public:
                       const amrex::MultiFab& Bz);
 #endif
 
+    void doRadiation(amrex::Real dt, amrex::Real cur_time);
+
+    void keepoldmomentum();
+    void keepoldmomentum_p(std::unique_ptr<WarpXParticleContainer>& pc);
+
+    void Dump_radiations();
     [[nodiscard]] int getSpeciesID (std::string product_str) const;
 
 protected:
@@ -480,6 +487,12 @@ protected:
     amrex::Real m_qed_schwinger_zmin = std::numeric_limits<amrex::Real>::lowest();
     amrex::Real m_qed_schwinger_zmax = std::numeric_limits<amrex::Real>::max();
 
+    /**
+     * For Radiation module
+     */
+    bool m_at_least_one_particle_radiate = 0;
+    int m_nspecies_radiate;
+
 #endif
 
 private:
@@ -527,6 +540,8 @@ private:
     void CheckQEDProductSpecies();
 #endif
 
-
+// Radiation Handler
+std::unique_ptr<RadiationHandler> m_p_radiation_handler;
+bool m_at_least_one_has_radiation = false;
 };
 #endif /*WARPX_ParticleContainer_H_*/

Source/Particles/MultiParticleContainer.cpp

@@ -122,7 +122,17 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core)
 
     // Setup particle collisions
     collisionhandler = std::make_unique<CollisionHandler>(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<RadiationHandler>();
+            break;
+        }
+    }
 }
 
 void
@@ -367,8 +377,7 @@ MultiParticleContainer::ReadParameters ()
                 m_qed_schwinger_threshold_poisson_gaussian);
             utils::parser::queryWithParser(
                 pp_qed_schwinger, "xmin", m_qed_schwinger_xmin);
-            utils::parser::queryWithParser(
-                pp_qed_schwinger, "xmax", m_qed_schwinger_xmax);
+
 #if defined(WARPX_DIM_3D)
             utils::parser::queryWithParser(
                 pp_qed_schwinger, "ymin", m_qed_schwinger_ymin);
@@ -384,7 +393,6 @@ MultiParticleContainer::ReadParameters ()
         initialized = true;
     }
 }
-
 WarpXParticleContainer&
 MultiParticleContainer::GetParticleContainerFromName (const std::string& name) const
 {
@@ -915,6 +923,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<WarpXParticleContainer>& 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];
+
+        });
+        }
+    }
+
 void
 MultiParticleContainer::doCollisions ( Real cur_time, amrex::Real dt )
 {
@@ -955,6 +1004,25 @@ void MultiParticleContainer::ScrapeParticles (const amrex::Vector<const amrex::M
 #endif
 }
 
+void MultiParticleContainer::doRadiation (const amrex::Real dt, amrex::Real cur_time)
+{
+    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);
+            }
+        }
+    }
+}
+
+
+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 ()
 {

Source/Particles/NamedComponentParticleContainer.H

@@ -177,6 +177,23 @@ public:
     /** Return the name-to-index map for the runtime-time integer components */
     [[nodiscard]] std::map<std::string, int> getParticleRuntimeiComps () const noexcept { return particle_runtime_icomps;}
 
+    /** Get the index of a real component
+     *
+     * @param name Name of the new component
+     */
+    int GetRealCompIndex (const std::string& name)
+    {
+        return particle_comps.at(name);
+    }
+     /** Initialize value at a component
+     *
+     * @param name Name of the new component
+     */
+    //    void InitializeComp (const std::string& name)
+    //{
+    //    particle_runtime_comps.at(name);
+    //}
+
 protected:
     std::map<std::string, int> particle_comps;
     std::map<std::string, int> particle_icomps;

Source/Particles/PhysicalParticleContainer.H

@@ -358,6 +358,8 @@ public:
     PairGenerationFilterFunc getPairGenerationFilterFunc ();
 #endif
 
+    bool has_radiation() const override {return m_do_radiation;}
+
 protected:
     std::string species_name;
     std::vector<std::unique_ptr<PlasmaInjector>> plasma_injectors;
@@ -386,6 +388,7 @@ protected:
     //When true PhysicalParticleContainer tries to use a pusher including
     //radiation reaction
     bool do_classical_radiation_reaction = false;
+    bool m_do_radiation = false;
 
     // A flag to enable saving of the previous timestep positions
     bool m_save_previous_position = false;

Source/Particles/PhysicalParticleContainer.cpp

@@ -347,6 +347,15 @@ 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==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
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(

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()

Source/Particles/Radiation/Make.package

@@ -0,0 +1 @@
+CEXE_sources += RadiationHandler.cpp

Source/Particles/Radiation/RadiationHandler.H

@@ -0,0 +1,86 @@
+
+/* 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/Pusher/GetAndSetPosition.H"
+
+
+#include <AMReX_REAL.H>
+#include <AMReX_Vector.H>
+#include <AMReX_Box.H>
+#include <AMReX_Vector.H>
+#include <AMReX_FArrayBox.H>
+#include <AMReX_Geometry.H>
+
+#include "WarpX.H"
+#include "Particles/Sorting/SortingUtils.H"
+#include <memory>
+#include <string>
+
+/* \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<WarpXParticleContainer>& pc, amrex::Real current_time);
+    void add_detector();
+    void gather_and_write_radiation(const std::string& filename);
+    void Integral_overtime();
+
+
+private:
+
+    // Frequency range
+    amrex::Vector<amrex::Real> m_omega_range;
+
+    // Dimensions of the detector
+    amrex::Vector<amrex::Real> m_theta_range;
+
+    int m_omega_points;
+
+    //put a detector
+    int m_get_a_detector;
+    // Resolution of the detector
+    amrex::Vector<amrex::Real>  m_d_theta;
+    double m_d_omega;
+    //The vector with the resolution of the detector
+    amrex::Vector<amrex::Real> m_d_det;
+    amrex::Real m_det_distance;
+
+    //Define the Fab with the datas
+    int ncomp;
+
+    amrex::Vector<int> m_det_pts;
+    amrex::Vector<int>  m_det_direction;
+    amrex::Vector<amrex::Real> m_pos_det;
+    amrex::Vector<amrex::Real> n;
+
+    //resolution of the grid of the detectir
+    amrex::Vector<amrex::Real> m_d_d;
+
+    //
+    amrex::Gpu::DeviceVector<amrex::Real> m_radiation_data;
+    amrex::Gpu::DeviceVector<amrex::Real> m_radiation_calculation;
+
+    //
+    amrex::Real dephas;
+
+    //
+    amrex::Vector<amrex::Vector<amrex::Real>> det_bornes;
+    amrex::Vector<amrex::Vector<amrex::Real>> det_pos;
+    amrex::Vector<amrex::Real> omega_calc;
+
+};
+#endif // WARPX_PARTICLES_RADIATION_H_