Merge branch 'dev' into update_perftest

1 files changed, 144 insertions, 47 deletions

diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp
index 1b644b543..9d39ec2f9 100644
--- a/Source/Particles/MultiParticleContainer.cpp
+++ b/Source/Particles/MultiParticleContainer.cpp
@@ -8,14 +8,11 @@
 
 using namespace amrex;
 
-constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr;
-
 MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core)
 {
     ReadParameters();
 
-    int n = WarpX::use_laser ? nspecies+1 : nspecies;
-    allcontainers.resize(n);
+    allcontainers.resize(nspecies + nlasers);
     for (int i = 0; i < nspecies; ++i) {
         if (species_types[i] == PCTypes::Physical) {
             allcontainers[i].reset(new PhysicalParticleContainer(amr_core, i, species_names[i]));
@@ -25,10 +22,46 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core)
         }
         allcontainers[i]->deposit_on_main_grid = deposit_on_main_grid[i];
     }
-    if (WarpX::use_laser) {
-	allcontainers[n-1].reset(new LaserParticleContainer(amr_core,n-1));
+    
+    for (int i = nspecies; i < nspecies+nlasers; ++i) {
+        allcontainers[i].reset(new LaserParticleContainer(amr_core,i, lasers_names[i-nspecies]));
+    }
+
+    pc_tmp.reset(new PhysicalParticleContainer(amr_core));
+
+    // Compute the number of species for which lab-frame data is dumped
+    // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer
+    // particle IDs in map_species_lab_diags.
+    map_species_boosted_frame_diags.resize(nspecies);
+    nspecies_boosted_frame_diags = 0;
+    for (int i=0; i<nspecies; i++){
+        auto& pc = allcontainers[i];
+        if (pc->do_boosted_frame_diags){
+            map_species_boosted_frame_diags[nspecies_boosted_frame_diags] = i;
+            do_boosted_frame_diags = 1;
+            nspecies_boosted_frame_diags += 1;
+        }
+    }
+
+    if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags)
+    {
+        for (int i = 0; i < nspecies_boosted_frame_diags; ++i)
+        {
+            int is = map_species_boosted_frame_diags[i];
+            allcontainers[is]->AddRealComp("xold");
+            allcontainers[is]->AddRealComp("yold");
+            allcontainers[is]->AddRealComp("zold");
+            allcontainers[is]->AddRealComp("uxold");
+            allcontainers[is]->AddRealComp("uyold");
+            allcontainers[is]->AddRealComp("uzold");
+        }
+        pc_tmp->AddRealComp("xold");
+        pc_tmp->AddRealComp("yold");
+        pc_tmp->AddRealComp("zold");
+        pc_tmp->AddRealComp("uxold");
+        pc_tmp->AddRealComp("uyold");
+        pc_tmp->AddRealComp("uzold");
     }
-    pc_tmp.reset(new PhysicalParticleContainer(amr_core));	
 }
 
 void
@@ -37,10 +70,10 @@ MultiParticleContainer::ReadParameters ()
     static bool initialized = false;
     if (!initialized)
     {
-	ParmParse pp("particles");
+        ParmParse pp("particles");
 
-	pp.query("nspecies", nspecies);
-	BL_ASSERT(nspecies >= 0);
+        pp.query("nspecies", nspecies);
+        BL_ASSERT(nspecies >= 0);
 
         if (nspecies > 0) {
             pp.getarr("species_names", species_names);
@@ -70,9 +103,19 @@ MultiParticleContainer::ReadParameters ()
                 }
             }
         }
-	pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr);
-	pp.query("l_lower_order_in_v", WarpX::l_lower_order_in_v);
-	initialized = true;
+
+        pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr);
+        pp.query("l_lower_order_in_v", WarpX::l_lower_order_in_v);
+
+        ParmParse ppl("lasers");
+        ppl.query("nlasers", nlasers);
+        BL_ASSERT(nlasers >= 0);
+        if (nlasers > 0) {
+            ppl.getarr("names", lasers_names);
+            BL_ASSERT(lasers_names.size() == nlasers);
+        }
+
+        initialized = true;
     }
 }
 
@@ -80,7 +123,7 @@ void
 MultiParticleContainer::AllocData ()
 {
     for (auto& pc : allcontainers) {
-	pc->AllocData();
+        pc->AllocData();
     }
     pc_tmp->AllocData();
 }
@@ -89,7 +132,7 @@ void
 MultiParticleContainer::InitData ()
 {
     for (auto& pc : allcontainers) {
-	pc->InitData();
+        pc->InitData();
     }
     pc_tmp->InitData();
 }
@@ -101,7 +144,7 @@ MultiParticleContainer::FieldGatherES (const Vector<std::array<std::unique_ptr<M
                                        const amrex::Vector<std::unique_ptr<amrex::FabArray<amrex::BaseFab<int> > > >& masks)
 {
     for (auto& pc : allcontainers) {
-	pc->FieldGatherES(E, masks);
+        pc->FieldGatherES(E, masks);
     }
 }
 
@@ -119,7 +162,7 @@ MultiParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<MultiF
     }
 
     for (auto& pc : allcontainers) {
-	pc->EvolveES(E, rho, t, dt);
+        pc->EvolveES(E, rho, t, dt);
     }
 
     for (unsigned i = 0; i < nlevs; i++) {
@@ -148,7 +191,7 @@ MultiParticleContainer::Evolve (int lev,
     if (rho) rho->setVal(0.0);
     for (auto& pc : allcontainers) {
 	pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz,
-                   rho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt);
+               rho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt);
     }    
 }
 
@@ -156,7 +199,7 @@ void
 MultiParticleContainer::PushXES (Real dt)
 {
     for (auto& pc : allcontainers) {
-	pc->PushXES(dt);
+        pc->PushXES(dt);
     }
 }
 
@@ -172,7 +215,7 @@ DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, bool local)
     }
 
     for (unsigned i = 0, n = allcontainers.size(); i < n; ++i) {
-	allcontainers[i]->DepositCharge(rho, true);
+        allcontainers[i]->DepositCharge(rho, true);
     }
 
     if (!local) {
@@ -201,7 +244,7 @@ MultiParticleContainer::FieldGather (int lev,
                                      const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
 {
     for (auto& pc : allcontainers) {
-	pc->FieldGather(lev, Ex, Ey, Ez, Bx, By, Bz);
+        pc->FieldGather(lev, Ex, Ey, Ez, Bx, By, Bz);
     }
 }
 
@@ -226,7 +269,7 @@ MultiParticleContainer::Evolve (int lev,
     if (crho) crho->setVal(0.0);
     for (auto& pc : allcontainers) {
 	pc->Evolve(lev, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz, cjx, cjy, cjz,
-                   rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt);
+               rho, crho, cEx, cEy, cEz, cBx, cBy, cBz, t, dt);
     }
 }
 
@@ -234,7 +277,7 @@ void
 MultiParticleContainer::PushX (Real dt)
 {
     for (auto& pc : allcontainers) {
-	pc->PushX(dt);
+        pc->PushX(dt);
     }
 }
 
@@ -244,7 +287,7 @@ MultiParticleContainer::PushP (int lev, Real dt,
                                const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz)
 {
     for (auto& pc : allcontainers) {
-       pc->PushP(lev, dt, Ex, Ey, Ez, Bx, By, Bz);
+        pc->PushP(lev, dt, Ex, Ey, Ez, Bx, By, Bz);
     }
 }
 
@@ -253,12 +296,12 @@ MultiParticleContainer::GetChargeDensity (int lev, bool local)
 {
     std::unique_ptr<MultiFab> rho = allcontainers[0]->GetChargeDensity(lev, true);
     for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) {
-	std::unique_ptr<MultiFab> rhoi = allcontainers[i]->GetChargeDensity(lev, true);
-	MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow());
+        std::unique_ptr<MultiFab> rhoi = allcontainers[i]->GetChargeDensity(lev, true);
+        MultiFab::Add(*rho, *rhoi, 0, 0, 1, rho->nGrow());
     }
     if (!local) {
-	const Geometry& gm = allcontainers[0]->Geom(lev);
-	rho->SumBoundary(gm.periodicity());
+        const Geometry& gm = allcontainers[0]->Geom(lev);
+        rho->SumBoundary(gm.periodicity());
     }
     return rho;
 }
@@ -267,7 +310,7 @@ void
 MultiParticleContainer::SortParticlesByCell ()
 {
     for (auto& pc : allcontainers) {
-	pc->SortParticlesByCell();
+        pc->SortParticlesByCell();
     }
 }
 
@@ -275,7 +318,7 @@ void
 MultiParticleContainer::Redistribute ()
 {
     for (auto& pc : allcontainers) {
-	pc->Redistribute();
+        pc->Redistribute();
     }
 }
 
@@ -283,7 +326,7 @@ void
 MultiParticleContainer::RedistributeLocal (const int num_ghost)
 {
     for (auto& pc : allcontainers) {
-	pc->Redistribute(0, 0, 0, num_ghost);
+        pc->Redistribute(0, 0, 0, num_ghost);
     }
 }
 
@@ -293,10 +336,10 @@ MultiParticleContainer::NumberOfParticlesInGrid(int lev) const
     const bool only_valid=true, only_local=true;
     Vector<long> r = allcontainers[0]->NumberOfParticlesInGrid(lev,only_valid,only_local);
     for (unsigned i = 1, n = allcontainers.size(); i < n; ++i) {
-	const auto& ri = allcontainers[i]->NumberOfParticlesInGrid(lev,only_valid,only_local);
-	for (unsigned j=0, m=ri.size(); j<m; ++j) {
-	    r[j] += ri[j];
-	}
+        const auto& ri = allcontainers[i]->NumberOfParticlesInGrid(lev,only_valid,only_local);
+        for (unsigned j=0, m=ri.size(); j<m; ++j) {
+            r[j] += ri[j];
+        }
     }
     ParallelDescriptor::ReduceLongSum(r.data(),r.size());
     return r;
@@ -306,7 +349,7 @@ void
 MultiParticleContainer::Increment (MultiFab& mf, int lev)
 {
     for (auto& pc : allcontainers) {
-	pc->Increment(mf,lev);
+        pc->Increment(mf,lev);
     }
 }
 
@@ -314,7 +357,7 @@ void
 MultiParticleContainer::SetParticleBoxArray (int lev, BoxArray& new_ba)
 {
     for (auto& pc : allcontainers) {
-	pc->SetParticleBoxArray(lev,new_ba);
+        pc->SetParticleBoxArray(lev,new_ba);
     }
 }
 
@@ -322,7 +365,7 @@ void
 MultiParticleContainer::SetParticleDistributionMap (int lev, DistributionMapping& new_dm)
 {
     for (auto& pc : allcontainers) {
-	pc->SetParticleDistributionMap(lev,new_dm);
+        pc->SetParticleDistributionMap(lev,new_dm);
     }
 }
 
@@ -330,7 +373,7 @@ void
 MultiParticleContainer::PostRestart ()
 {
     for (auto& pc : allcontainers) {
-	pc->PostRestart();
+        pc->PostRestart();
     }
     pc_tmp->PostRestart();
 }
@@ -346,16 +389,25 @@ MultiParticleContainer
 
     BL_PROFILE("MultiParticleContainer::GetLabFrameData");
     
-    for (int i = 0; i < nspecies; ++i)
-    {
-        WarpXParticleContainer* pc = allcontainers[i].get();
+    // Loop over particle species
+    for (int i = 0; i < nspecies_boosted_frame_diags; ++i){
+        int isp = map_species_boosted_frame_diags[i];
+        WarpXParticleContainer* pc = allcontainers[isp].get();
         WarpXParticleContainer::DiagnosticParticles diagnostic_particles;
         pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles);
-
-        for (int lev = 0; lev <= pc->finestLevel(); ++lev)
-        {
-            for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it)
-            {
+        // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData
+        // where "lev" is the AMR level and "index" is a [grid index][tile index] pair.
+        
+        // Loop over AMR levels
+        for (int lev = 0; lev <= pc->finestLevel(); ++lev){
+            // Loop over [grid index][tile index] pairs 
+            // and Fills parts[species number i] with particle data from all grids and 
+            // tiles in diagnostic_particles. parts contains particles from all 
+            // AMR levels indistinctly.
+            for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){
+                // it->first is the [grid index][tile index] key
+                // it->second is the corresponding 
+                // WarpXParticleContainer::DiagnosticParticleData value
                 parts[i].GetRealData(DiagIdx::w).insert(  parts[i].GetRealData(DiagIdx::w  ).end(),
                                                           it->second.GetRealData(DiagIdx::w  ).begin(),
                                                           it->second.GetRealData(DiagIdx::w  ).end());
@@ -387,3 +439,48 @@ MultiParticleContainer
         }
     }
 }
+
+/* \brief Continuous injection for particles initially outside of the domain.
+ * \param injection_box: Domain where new particles should be injected.
+ * Loop over all WarpXParticleContainer in MultiParticleContainer and 
+ * calls virtual function ContinuousInjection.
+ */
+void
+MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const
+{
+    for (int i=0; i<nspecies+nlasers; i++){
+        auto& pc = allcontainers[i];
+        if (pc->do_continuous_injection){
+            pc->ContinuousInjection(injection_box);
+        }
+    }
+}
+
+/* \brief Update position of continuous injection parameters.
+ * \param dt: simulation time step (level 0)
+ * All classes inherited from WarpXParticleContainer do not have 
+ * a position to update (PhysicalParticleContainer does not do anything).
+ */
+void
+MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const
+{
+    for (int i=0; i<nspecies+nlasers; i++){
+        auto& pc = allcontainers[i];
+        if (pc->do_continuous_injection){
+            pc->UpdateContinuousInjectionPosition(dt);
+        }
+    }
+}
+
+int
+MultiParticleContainer::doContinuousInjection() const
+{
+    int warpx_do_continuous_injection = 0;
+    for (int i=0; i<nspecies+nlasers; i++){
+        auto& pc = allcontainers[i];
+        if (pc->do_continuous_injection){
+            warpx_do_continuous_injection = 1;
+        }
+    }
+    return warpx_do_continuous_injection;
+}