Obtain exact time (real) when particles hit boundaries + test (#4695)

* adding the time

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update particle_containers.py

* Update ParticleBoundaryBuffer.cpp

* Update PICMI_inputs_rz.py

* Update PICMI_inputs_rz.py

* Update particle_containers.py

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* fix issue

* fix issue

* issue fixed

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update analysis_rz.py

* Update analysis_rz_filter.py

* Update analysis.py

* Update PICMI_inputs_scrape.py

* Update PICMI_inputs_reflection.py

* Update particle_containers.py

* Update PICMI_inputs_reflection.py

* Update PICMI_inputs_scrape.py

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update ParticleBoundaryBuffer.cpp

* Update Source/Particles/ParticleBoundaryBuffer.cpp

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Update Source/Particles/ParticleBoundaryBuffer.cpp

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Update ParticleBoundaryBuffer.cpp

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

Author: 3 people

Docs/source/usage/parameters.rst

@@ -2752,8 +2752,10 @@ The data collected at each boundary is written out to a subdirectory of the diag
 By default, all of the collected particle data is written out at the end of the simulation. Optionally, the ``<diag_name>.intervals`` parameter can be given to specify writing out the data more often.
 This can be important if a large number of particles are lost, avoiding filling up memory with the accumulated lost particle data.
 
-In addition to their usual attributes, the saved particles have an integer attribute ``timestamp``, which
-indicates the PIC iteration at which each particle was absorbed at the boundary.
+In addition to their usual attributes, the saved particles have an integer attribute ``step_scraped``, which
+indicates the PIC iteration at which each particle was absorbed at the boundary,
+and a real attribute ``time_scraped``, which indicates the exact time calculated when each
+particle hits the EB.
 
 ``BoundaryScrapingDiagnostics`` can be used with ``<diag_name>.<species>.random_fraction``, ``<diag_name>.<species>.uniform_stride``, and ``<diag_name>.<species>.plot_filter_function``, which have the same behavior as for ``FullDiagnostics``. For ``BoundaryScrapingDiagnostics``, these filters are applied at the time the data is written to file. An implication of this is that more particles may initially be accumulated in memory than are ultimately written. ``t`` in ``plot_filter_function`` refers to the time the diagnostic is written rather than the time the particle crossed the boundary.
 

Examples/Tests/point_of_contact_EB/analysis.py

@@ -26,26 +26,35 @@
 ts_scraping = OpenPMDTimeSeries('./diags/diag2/particles_at_eb/')
 
 it=ts_scraping.iterations
-x,y,z=ts_scraping.get_particle( ['x','y','z'], species='electron', iteration=it )
+step_scraped, time_scraped, x, y, z=ts_scraping.get_particle( ['stepScraped','timeScraped','x','y','z'], species='electron', iteration=it )
+time_scraped_reduced=time_scraped[0]*1e10
 
 # Analytical results calculated
 x_analytic=-0.1983
 y_analytic=0.02584
 z_analytic=0.0000
 
+#result obtained by analysis of simulations
+step=3
+time_reduced=3.58
+
 print('NUMERICAL coordinates of the point of contact:')
-print('x=%5.4f, y=%5.4f, z=%5.4f' % (x[0], y[0], z[0]))
+print('step_scraped=%d, time_stamp=%5.4f e-10, x=%5.4f, y=%5.4f, z=%5.4f' % (step_scraped[0],time_reduced,x[0], y[0], z[0]))
 print('\n')
 print('ANALYTICAL coordinates of the point of contact:')
-print('x=%5.4f, y=%5.4f, z=%5.4f' % (x_analytic, y_analytic, z_analytic))
+print('step_scraped=%d, time_stamp=%5.4f e-10, x=%5.4f, y=%5.4f, z=%5.4f' % (step, time_reduced, x_analytic, y_analytic, z_analytic))
 
 tolerance=0.001
+tolerance_t=0.003
 print("tolerance = "+ str(tolerance *100) + '%')
+print("tolerance for the time = "+ str(tolerance_t *100) + '%')
 
+diff_step=np.abs((step_scraped[0]-step)/step)
+diff_time=np.abs((time_scraped_reduced-time_reduced)/time_reduced)
 diff_x=np.abs((x[0]-x_analytic)/x_analytic)
 diff_y=np.abs((y[0]-y_analytic)/y_analytic)
 
 print("percentage error for x = %5.4f %%" %(diff_x *100))
 print("percentage error for y = %5.4f %%" %(diff_y *100))
 
-assert (diff_x < tolerance) and (diff_y < tolerance) and (np.abs(z[0]) < 1e-8) , 'Test point_of_contact did not pass'
+assert (diff_x < tolerance) and (diff_y < tolerance) and (np.abs(z[0]) < 1e-8) and (diff_step < 1e-8) and (diff_time < tolerance_t) , 'Test point_of_contact did not pass'

Examples/Tests/scraping/analysis_rz.py

@@ -53,8 +53,8 @@ def n_remaining_particles( iteration ):
     w, = ts_full.get_particle(['w'], iteration=iteration)
     return len(w)
 def n_scraped_particles( iteration ):
-    timestamp = ts_scraping.get_particle( ['timestamp'], iteration=ts_scraping.iterations[0] )
-    return (timestamp <= iteration).sum()
+    step_scraped = ts_scraping.get_particle( ['stepScraped'], iteration=ts_scraping.iterations[0] )
+    return (step_scraped <= iteration).sum()
 n_remaining = np.array([ n_remaining_particles(iteration) for iteration in ts_full.iterations ])
 n_scraped = np.array([ n_scraped_particles(iteration) for iteration in ts_full.iterations ])
 n_total = n_remaining[0]

Examples/Tests/scraping/analysis_rz_filter.py

@@ -50,8 +50,8 @@ def n_remaining_particles( iteration ):
     w, = ts_full.get_particle(['w'], iteration=iteration)
     return len(w)
 def n_scraped_particles( iteration ):
-    timestamp = ts_scraping.get_particle( ['timestamp'], iteration=ts_scraping.iterations[0] )
-    return (timestamp <= iteration).sum()
+    step_scraped = ts_scraping.get_particle( ['stepScraped'], iteration=ts_scraping.iterations[0] )
+    return (step_scraped <= iteration).sum()
 def n_scraped_z_leq_zero( iteration ):
     z_pos, = ts_scraping.get_particle( ['z'], iteration=ts_scraping.iterations[0] )
     return (z_pos <= 0).sum()

Python/pywarpx/particle_containers.py

@@ -772,9 +772,9 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level)
                 form x/y/z_hi/lo or eb.
 
             comp_name      : str
-                The component of the array data that will be returned. If
-                "step_scraped" the special attribute holding the timestep at
-                which a particle was scraped will be returned.
+                The component of the array data that will be returned.
+                "x", "y", "z", "ux", "uy", "uz", "w"
+                "step_scraped","time_scraped", "nx", "ny", "nz"
 
             level          : int
                 Which AMR level to retrieve scraped particle data from.
@@ -785,18 +785,20 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level)
             species_name, self._get_boundary_number(boundary)
         )
         data_array = []
-        if comp_name == 'step_scraped':
-            # the step scraped is always the final integer component
-            comp_idx = part_container.num_int_comps - 1
+        #loop over the real attributes
+        if comp_name in part_container.real_comp_names:
+            comp_idx = part_container.real_comp_names[comp_name]
             for ii, pti in enumerate(libwarpx.libwarpx_so.BoundaryBufferParIter(part_container, level)):
+                soa = pti.soa()
+                data_array.append(xp.array(soa.get_real_data(comp_idx), copy=False))
+        #loop over the integer attributes
+        elif comp_name in part_container.int_comp_names:
+             comp_idx = part_container.int_comp_names[comp_name]
+             for ii, pti in enumerate(libwarpx.libwarpx_so.BoundaryBufferParIter(part_container, level)):
                 soa = pti.soa()
                 data_array.append(xp.array(soa.get_int_data(comp_idx), copy=False))
         else:
-            container_wrapper = ParticleContainerWrapper(species_name)
-            comp_idx = container_wrapper.particle_container.get_comp_index(comp_name)
-            for ii, pti in enumerate(libwarpx.libwarpx_so.BoundaryBufferParIter(part_container, level)):
-                soa = pti.soa()
-                data_array.append(xp.array(soa.get_real_data(comp_idx), copy=False))
+            raise RuntimeError('Name %s not found' %comp_name)
         return data_array
 
 

Source/Diagnostics/BTDiagnostics.cpp

@@ -1299,7 +1299,7 @@ BTDiagnostics::InterleaveBufferAndSnapshotHeader ( std::string buffer_Header_pat
     BTDPlotfileHeaderImpl buffer_HeaderImpl(buffer_Header_path);
     buffer_HeaderImpl.ReadHeaderData();
 
-    // Update timestamp of snapshot with recently flushed buffer
+    // Update step_scraped of snapshot with recently flushed buffer
     snapshot_HeaderImpl.set_time( buffer_HeaderImpl.time() );
     snapshot_HeaderImpl.set_timestep( buffer_HeaderImpl.timestep() );
 

Source/Particles/ParticleBoundaryBuffer.cpp

@@ -45,7 +45,8 @@ struct IsOutsideDomainBoundary {
 };
 
 struct FindEmbeddedBoundaryIntersection {
-    const int m_index;
+    const int m_step_index;
+    const int m_time_index;
     const int m_step;
     const amrex::Real m_dt;
     amrex::Array4<const amrex::Real> m_phiarr;
@@ -69,9 +70,6 @@ struct FindEmbeddedBoundaryIntersection {
             dst.m_runtime_idata[j][dst_i] = src.m_runtime_idata[j][src_i];
         }
 
-        // Also record the integer timestep on the destination
-        dst.m_runtime_idata[m_index][dst_i] = m_step;
-
         // Modify the position of the destination particle:
         // Move it to the point of intersection with the embedded boundary
         // (which is found by using a bisection algorithm)
@@ -90,7 +88,6 @@ struct FindEmbeddedBoundaryIntersection {
         amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> plo = m_plo;
 
         // Bisection algorithm to find the point where phi(x,y,z)=0 (i.e. on the embedded boundary)
-
         amrex::Real dt_fraction = amrex::bisect( 0.0, 1.0,
             [=] (amrex::Real dt_frac) {
                 int i, j, k;
@@ -102,6 +99,10 @@ struct FindEmbeddedBoundaryIntersection {
                 return phi_value;
             } );
 
+        // Also record the real time on the destination
+        dst.m_runtime_idata[m_step_index][dst_i] = m_step;
+        dst.m_runtime_rdata[m_time_index][dst_i] = m_step*m_dt + (1- dt_fraction)*m_dt;
+
         // Now that dt_fraction has be obtained (with bisect)
         // Save the corresponding position of the particle at the boundary
         amrex::ParticleReal x_temp=xp, y_temp=yp, z_temp=zp;
@@ -129,8 +130,10 @@ struct FindEmbeddedBoundaryIntersection {
 };
 
 struct CopyAndTimestamp {
-    int m_index;
+    int m_step_index;
+    int m_time_index;
     int m_step;
+    const amrex::Real m_dt;
 
     template <typename DstData, typename SrcData>
     AMREX_GPU_HOST_DEVICE
@@ -147,10 +150,12 @@ struct CopyAndTimestamp {
         for (int j = 0; j < src.m_num_runtime_int; ++j) {
             dst.m_runtime_idata[j][dst_i] = src.m_runtime_idata[j][src_i];
         }
-        dst.m_runtime_idata[m_index][dst_i] = m_step;
+        dst.m_runtime_idata[m_step_index][dst_i] = m_step;
+        dst.m_runtime_rdata[m_time_index][dst_i] = m_step*m_dt;
     }
 };
 
+
 ParticleBoundaryBuffer::ParticleBoundaryBuffer ()
 {
     m_particle_containers.resize(numBoundaries());
@@ -328,7 +333,8 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc,
                 if (!buffer[i].isDefined())
                 {
                     buffer[i] = pc.make_alike<amrex::PinnedArenaAllocator>();
-                    buffer[i].AddIntComp("timestamp", false);
+                    buffer[i].AddIntComp("step_scraped", false);
+                    buffer[i].AddRealComp("time_scraped", false);
                 }
                 auto& species_buffer = buffer[i];
                 for (int lev = 0; lev < pc.numLevels(); ++lev)
@@ -368,12 +374,17 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc,
                         }
                         {
                           WARPX_PROFILE("ParticleBoundaryBuffer::gatherParticles::filterAndTransform");
-                          const int timestamp_index = ptile_buffer.NumRuntimeIntComps()-1;
-                          const int timestep = warpx_instance.getistep(0);
+                          auto& warpx = WarpX::GetInstance();
+                          const auto dt = warpx.getdt(pti.GetLevel());
+                          auto string_to_index_intcomp = buffer[i].getParticleRuntimeiComps();
+                          const int step_scraped_index = string_to_index_intcomp.at("step_scraped");
+                          auto string_to_index_realcomp = buffer[i].getParticleRuntimeComps();
+                          const int time_scraped_index = string_to_index_realcomp.at("time_scraped");
+                          const int step = warpx_instance.getistep(0);
 
                           amrex::filterAndTransformParticles(ptile_buffer, ptile,
                                                              predicate,
-                                                             CopyAndTimestamp{timestamp_index, timestep},
+                                                             CopyAndTimestamp{step_scraped_index,time_scraped_index, step,dt},
                                                              0, dst_index);
                         }
                     }
@@ -393,7 +404,8 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc,
         if (!buffer[i].isDefined())
         {
             buffer[i] = pc.make_alike<amrex::PinnedArenaAllocator>();
-            buffer[i].AddIntComp("timestamp", false);
+            buffer[i].AddIntComp("step_scraped", false);
+            buffer[i].AddRealComp("time_scraped", false);
         }
         auto& species_buffer = buffer[i];
         for (int lev = 0; lev < pc.numLevels(); ++lev)
@@ -444,16 +456,18 @@ void ParticleBoundaryBuffer::gatherParticles (MultiParticleContainer& mypc,
                   WARPX_PROFILE("ParticleBoundaryBuffer::gatherParticles::resize_eb");
                   ptile_buffer.resize(dst_index + amrex::get<0>(reduce_data.value()));
                 }
-
-                const int timestamp_index = ptile_buffer.NumRuntimeIntComps()-1;
-                const int timestep = warpx_instance.getistep(0);
                 auto& warpx = WarpX::GetInstance();
                 const auto dt = warpx.getdt(pti.GetLevel());
+                auto string_to_index_intcomp = buffer[i].getParticleRuntimeiComps();
+                const int step_scraped_index = string_to_index_intcomp.at("step_scraped");
+                auto string_to_index_realcomp = buffer[i].getParticleRuntimeComps();
+                const int time_scraped_index = string_to_index_realcomp.at("time_scraped");
+                const int step = warpx_instance.getistep(0);
 
                 {
                   WARPX_PROFILE("ParticleBoundaryBuffer::gatherParticles::filterTransformEB");
                   amrex::filterAndTransformParticles(ptile_buffer, ptile, predicate,
-                                                     FindEmbeddedBoundaryIntersection{timestamp_index, timestep, dt, phiarr, dxi, plo}, 0, dst_index);
+                                                     FindEmbeddedBoundaryIntersection{step_scraped_index,time_scraped_index, step, dt, phiarr, dxi, plo}, 0, dst_index);
                 }
             }
         }

Source/Python/Particles/PinnedMemoryParticleContainer.cpp

@@ -15,4 +15,17 @@ void init_PinnedMemoryParticleContainer (py::module& m)
         PinnedMemoryParticleContainer,
         amrex::ParticleContainerPureSoA<PIdx::nattribs, 0, amrex::PinnedArenaAllocator>
     > pmpc (m, "PinnedMemoryParticleContainer");
+    pmpc
+        .def_property_readonly("real_comp_names",
+            [](PinnedMemoryParticleContainer& pc)
+            {
+                return pc.getParticleComps();
+            }
+        )
+        .def_property_readonly("int_comp_names",
+            [](PinnedMemoryParticleContainer& pc)
+            {
+                return pc.getParticleiComps();
+            }
+        );
 }