Now physical particles have Tau initialization via QS engine

1 files changed, 92 insertions, 60 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index da4e3afc2..ae226015d 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
     pp.query("do_splitting", do_splitting);
     pp.query("split_type", split_type);
     pp.query("do_continuous_injection", do_continuous_injection);
-    // Whether to plot back-transformed (lab-frame) diagnostics 
+    // Whether to plot back-transformed (lab-frame) diagnostics
     // for this species.
     pp.query("do_boosted_frame_diags", do_boosted_frame_diags);
 
@@ -50,9 +50,14 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
 #ifdef WARPX_QED
     //Add real component if QED is enabled
     pp.query("do_qed", do_qed);
-    if(do_qed){
+    if(do_qed)
         AddRealComp("tau");
-    }
+
+    //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled
+    if(do_qed)
+        pp.query("do_qed_quantum_sync", do_qed_quantum_sync);
+
+    //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!!
 #endif
 
     //variable to set plot_flags size
@@ -63,7 +68,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
 #ifdef WARPX_QED
     if(do_qed){
         // plot_flag will have an entry for the optical depth
-        plot_flag_size ++;
+        plot_flag_size++;
     }
 #endif
     //_______________________________
@@ -104,7 +109,14 @@ void PhysicalParticleContainer::InitData()
     // Init ionization module here instead of in the PhysicalParticleContainer
     // constructor because dt is required
     if (do_field_ionization) {InitIonizationModule();}
+
     AddParticles(0); // Note - add on level 0
+
+#ifdef WARPX_QED
+    if(do_qed_quantum_sync)
+        InitTauQuantumSync();
+#endif
+
     Redistribute();  // We then redistribute
 }
 
@@ -157,7 +169,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real
 void
 PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
                                            Real x_rms, Real y_rms, Real z_rms,
-                                           Real q_tot, long npart, 
+                                           Real q_tot, long npart,
                                            int do_symmetrize) {
 
     const Geometry& geom     = m_gdb->Geom(0);
@@ -169,7 +181,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
     std::normal_distribution<double> distz(z_m, z_rms);
 
     if (ParallelDescriptor::IOProcessor()) {
-        // If do_symmetrize, create 4x fewer particles, and 
+        // If do_symmetrize, create 4x fewer particles, and
         // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y)
         if (do_symmetrize){
             npart /= 4;
@@ -402,11 +414,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
         const int tile_id = mfi.LocalTileIndex();
 
         // Max number of new particles, if particles are created in the whole
-        // overlap_box. All of them are created, and invalid ones are then 
+        // overlap_box. All of them are created, and invalid ones are then
         // discaded
         int max_new_particles = overlap_box.numPts() * num_ppc;
 
-        // If refine injection, build pointer dp_cellid that holds pointer to 
+        // If refine injection, build pointer dp_cellid that holds pointer to
         // array of refined cell IDs.
         Vector<int> cellid_v;
         if (refine_injection and lev == 0)
@@ -461,7 +473,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
         if (do_field_ionization) {
             pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size;
         }
-        
+
         const GpuArray<Real,AMREX_SPACEDIM> overlap_corner
             {AMREX_D_DECL(overlap_realbox.lo(0),
                           overlap_realbox.lo(1),
@@ -473,9 +485,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
 	bool loc_do_field_ionization = do_field_ionization;
 	int loc_ionization_initial_level = ionization_initial_level;
 
-        // Loop over all new particles and inject them (creates too many 
+        // Loop over all new particles and inject them (creates too many
         // particles, in particular does not consider xmin, xmax etc.).
-        // The invalid ones are given negative ID and are deleted during the 
+        // The invalid ones are given negative ID and are deleted during the
         // next redistribute.
         amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept
         {
@@ -642,7 +654,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox)
             p.pos(1) = z;
 #endif
         }, shared_mem_bytes);
-    			 
+
         if (cost) {
             wt = (amrex::second() - wt) / tile_box.d_numPts();
             Array4<Real> const& costarr = cost->array(mfi);
@@ -880,7 +892,7 @@ PhysicalParticleContainer::FieldGather (int lev,
     MultiFab* cost = WarpX::getCosts(lev);
 
 #ifdef _OPENMP
-#pragma omp parallel 
+#pragma omp parallel
 #endif
     {
 #ifdef _OPENMP
@@ -930,7 +942,7 @@ PhysicalParticleContainer::FieldGather (int lev,
             //
             int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
             FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, 
+                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
                         Ex.nGrow(), e_is_nodal,
                         0, np, thread_num, lev, lev);
 
@@ -964,11 +976,11 @@ PhysicalParticleContainer::Evolve (int lev,
     BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg);
     BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp);
     BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition);
-    
+
     const std::array<Real,3>& dx = WarpX::CellSize(lev);
     const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0));
 
-    // Get instances of NCI Godfrey filters 
+    // Get instances of NCI Godfrey filters
     const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz;
     const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez;
 
@@ -993,9 +1005,9 @@ PhysicalParticleContainer::Evolve (int lev,
                 tmp_particle_data[lev][index][i].resize(np);
         }
     }
-    
+
 #ifdef _OPENMP
-#pragma omp parallel 
+#pragma omp parallel
 #endif
     {
 #ifdef _OPENMP
@@ -1194,7 +1206,7 @@ PhysicalParticleContainer::Evolve (int lev,
             }
 
             const long np_current = (cjx) ? nfine_current : np;
-            
+
             //
             // copy data from particle container to temp arrays
             //
@@ -1210,14 +1222,14 @@ PhysicalParticleContainer::Evolve (int lev,
                 } else {
                     ion_lev = nullptr;
                 }
-                DepositCharge(pti, wp, ion_lev, rho, 0, 0, 
+                DepositCharge(pti, wp, ion_lev, rho, 0, 0,
                               np_current, thread_num, lev, lev);
                 if (has_buffer){
                     DepositCharge(pti, wp, ion_lev, crho, 0, np_current,
                                   np-np_current, thread_num, lev, lev-1);
                 }
             }
-            
+
             if (! do_not_push)
             {
                 const long np_gather = (cEx) ? nfine_gather : np;
@@ -1229,7 +1241,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 //
                 BL_PROFILE_VAR_START(blp_fg);
                 FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                            exfab, eyfab, ezfab, bxfab, byfab, bzfab, 
+                            exfab, eyfab, ezfab, bxfab, byfab, bzfab,
                             Ex.nGrow(), e_is_nodal,
                             0, np_gather, thread_num, lev, lev);
 
@@ -1245,7 +1257,7 @@ PhysicalParticleContainer::Evolve (int lev,
                     FArrayBox const* cbxfab = &(*cBx)[pti];
                     FArrayBox const* cbyfab = &(*cBy)[pti];
                     FArrayBox const* cbzfab = &(*cBz)[pti];
-                    
+
                     if (WarpX::use_fdtd_nci_corr)
                     {
 #if (AMREX_SPACEDIM == 2)
@@ -1283,13 +1295,13 @@ PhysicalParticleContainer::Evolve (int lev,
                         eyeli = filtered_Ey.elixir();
                         nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box());
                         ceyfab = &filtered_Ey;
-                        
+
                         // Filter Bx
                         filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag));
                         bxeli = filtered_Bx.elixir();
                         nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box());
                         cbxfab = &filtered_Bx;
-                        
+
                         // Filter Bz
                         filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag));
                         bzeli = filtered_Bz.elixir();
@@ -1297,14 +1309,14 @@ PhysicalParticleContainer::Evolve (int lev,
                         cbzfab = &filtered_Bz;
 #endif
                     }
-                    
+
                     // Field gather for particles in gather buffers
                     e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal();
-                    FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, 
+                    FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
                                 cexfab, ceyfab, cezfab,
                                 cbxfab, cbyfab, cbzfab,
-                                cEx->nGrow(), e_is_nodal, 
-                                nfine_gather, np-nfine_gather, 
+                                cEx->nGrow(), e_is_nodal,
+                                nfine_gather, np-nfine_gather,
                                 thread_num, lev, lev-1);
                 }
 
@@ -1314,7 +1326,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 // Particle Push
                 //
                 BL_PROFILE_VAR_START(blp_ppc_pp);
-                PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], 
+                PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num],
                        m_giv[thread_num], dt);
                 BL_PROFILE_VAR_STOP(blp_ppc_pp);
 
@@ -1328,7 +1340,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 } else {
                     ion_lev = nullptr;
                 }
-                
+
                 // Deposit inside domains
                 DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz,
                                0, np_current, thread_num,
@@ -1348,7 +1360,7 @@ PhysicalParticleContainer::Evolve (int lev,
                 pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]);
                 BL_PROFILE_VAR_STOP(blp_copy);
             }
-            
+
             if (rho) {
                 // Deposit charge after particle push, in component 1 of MultiFab rho.
                 int* AMREX_RESTRICT ion_lev;
@@ -1417,7 +1429,7 @@ PhysicalParticleContainer::SplitParticles(int lev)
         for(int i=0; i<np; i++){
             auto& p = particles[i];
             if (p.id() == DoSplitParticleID){
-                // If particle is tagged, split it and put the 
+                // If particle is tagged, split it and put the
                 // split particles in local arrays psplit_x etc.
                 np_split_to_add += np_split;
 #if (AMREX_SPACEDIM==2)
@@ -1514,11 +1526,11 @@ PhysicalParticleContainer::SplitParticles(int lev)
     }
 	// Add local arrays psplit_x etc. to the temporary
 	// particle container pctmp_split. Split particles
-	// are tagged with p.id()=NoSplitParticleID so that 
+	// are tagged with p.id()=NoSplitParticleID so that
 	// they are not re-split when entering a higher level
 	// AddNParticles calls Redistribute, so that particles
 	// in pctmp_split are in the proper grids and tiles
-	pctmp_split.AddNParticles(lev, 
+	pctmp_split.AddNParticles(lev,
                               np_split_to_add,
                               psplit_x.dataPtr(),
                               psplit_y.dataPtr(),
@@ -1570,12 +1582,12 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
     if (do_field_ionization){
         ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
     }
-    
+
     // Loop over the particles and update their momentum
     const Real q = this->charge;
     const Real m = this-> mass;
     if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
-        amrex::ParallelFor( 
+        amrex::ParallelFor(
             pti.numParticles(),
             [=] AMREX_GPU_DEVICE (long i) {
                 Real qp = q;
@@ -1624,7 +1636,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
         int thread_num = omp_get_thread_num();
 #else
         int thread_num = 0;
-#endif      
+#endif
         for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti)
         {
             const Box& box = pti.validbox();
@@ -1664,7 +1676,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
 
             int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal();
             FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp,
-                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, 
+                        &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab,
                         Ex.nGrow(), e_is_nodal,
                         0, np, thread_num, lev, lev);
 
@@ -1712,7 +1724,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp,
     Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr();
     Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr();
     Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr();
-    
+
     const auto np = pti.numParticles();
     const auto lev = pti.GetLevel();
     const auto index = pti.GetPairIndex();
@@ -1728,7 +1740,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp,
                      xpold[i]=xp[i];
                      ypold[i]=yp[i];
                      zpold[i]=zp[i];
-            
+
                      uxpold[i]=uxp[i];
                      uypold[i]=uyp[i];
                      uzpold[i]=uzp[i];
@@ -1768,7 +1780,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
     slice_box.setHi(direction, z_max);
 
     diagnostic_particles.resize(finestLevel()+1);
-    
+
     for (int lev = 0; lev < nlevs; ++lev) {
 
         const Real* dx  = Geom(lev).CellSize();
@@ -1815,7 +1827,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
                 auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold];
 
                 const long np = pti.numParticles();
-                
+
                 Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c;
                 Real inv_c2 = 1.0/PhysConst::c/PhysConst::c;
 
@@ -1850,7 +1862,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real
                     Real uzp = uz_old_p  *weight_old + uz_new_p  *weight_new;
 
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]);
-                    
+
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp);
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp);
                     diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp);
@@ -1875,7 +1887,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box)
     AddPlasma(lev, injection_box);
 }
 
-/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, 
+/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab,
  * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp.
  * \param Exp-Bzp: fields on particles.
  * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti
@@ -1885,7 +1897,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box)
  * \param np_to_gather: number of particles onto which fields are gathered
  * \param thread_num: if using OpenMP, thread number
  * \param lev: level on which particles are located
- * \param gather_lev: level from which particles gather fields (lev-1) for 
+ * \param gather_lev: level from which particles gather fields (lev-1) for
           particles in buffers.
  */
 void
@@ -1912,14 +1924,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) ||
                                      (gather_lev==(lev  )),
                                      "Gather buffers only work for lev-1");
-    
+
     // If no particles, do not do anything
     if (np_to_gather == 0) return;
     // Get cell size on gather_lev
     const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0));
     // Set staggering shift depending on e_is_nodal
     const Real stagger_shift = e_is_nodal ? 0.0 : 0.5;
-    
+
     // Get box from which field is gathered.
     // If not gathering from the finest level, the box is coarsened.
     Box box;
@@ -1929,26 +1941,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti,
         const IntVect& ref_ratio = WarpX::RefRatio(gather_lev);
         box = amrex::coarsen(pti.tilebox(),ref_ratio);
     }
-    
+
     // Add guard cells to the box.
     box.grow(ngE);
-    
+
     const Array4<const Real>& ex_arr = exfab->array();
     const Array4<const Real>& ey_arr = eyfab->array();
     const Array4<const Real>& ez_arr = ezfab->array();
     const Array4<const Real>& bx_arr = bxfab->array();
     const Array4<const Real>& by_arr = byfab->array();
     const Array4<const Real>& bz_arr = bzfab->array();
-    
+
     const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset;
     const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset;
     const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset;
-    
+
     // Lower corner of tile box physical domain
     const std::array<Real, 3>& xyzmin = WarpX::LowerCorner(box, gather_lev);
-    
+
     const Dim3 lo = lbound(box);
-    
+
     // Depending on l_lower_in_v and WarpX::nox, call
     // different versions of template function doGatherShapeN
     if (WarpX::l_lower_order_in_v){
@@ -2031,10 +2043,10 @@ void PhysicalParticleContainer::InitIonizationModule ()
     }
     // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2))
     // For now, we assume l=0 and m=0.
-    // The approximate expressions are used, 
+    // The approximate expressions are used,
     // without Gamma function
     Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e;
-    Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * 
+    Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e *
         std::pow(PhysConst::alpha,4)/PhysConst::r_e;
     Real UH = table_ionization_energies[0];
     Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.;
@@ -2049,18 +2061,18 @@ void PhysicalParticleContainer::InitIonizationModule ()
         Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff));
         adk_power[i] = -(2*n_eff - 1);
         Real Uion = ionization_energies[i];
-        adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) 
+        adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) )
             * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1);
         adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea;
     }
 }
 
 /* \brief create mask of ionized particles (1 if ionized, 0 otherwise)
- * 
+ *
  * \param mfi: tile or grid
  * \param lev: MR level
  * \param ionization_mask: Array with as many elements as particles in mfi.
- * This function initialized the array, and set each element to 1 or 0 
+ * This function initialized the array, and set each element to 1 or 0
  * depending on whether the particle is ionized or not.
  */
 void
@@ -2106,7 +2118,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
 
     // Loop over all particles in grid/tile. If ionized, set mask to 1
     // and increment ionization level.
-    ParallelFor( 
+    ParallelFor(
         np,
         [=] AMREX_GPU_DEVICE (long i) {
             // Get index of ionization_level
@@ -2138,3 +2150,23 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
         }
     );
 }
+
+#ifdef WARPX_QED
+// A function to initialize the Tau component according to the QS engine
+void PhysicalParticleContainer::InitTauQuantumSync()
+{
+    BL_PROFILE("PhysicalParticleContainer::InitTauQuantumSync");
+//Looping over all the particles
+    int num_levels = finestLevel() + 1;
+    for (int lev=0; lev < num_levels; ++lev)
+        for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){
+            auto taus = pti.GetAttribs(particle_comps["tau"]).dataPtr();
+            amrex::ParallelFor(
+                pti.numParticles(),
+                [=] AMREX_GPU_DEVICE (long i) {
+                    taus[i] = warpx_quantum_sync_engine::get_optical_depth();
+                }
+                );
+    }
+}
+#endif