Merge branch 'development' into dt_fraction

.azure-pipelines.yml

@@ -41,7 +41,7 @@ jobs:
         WARPX_CI_EB: 'TRUE'
 
   # default: 60; maximum: 360
-  timeoutInMinutes: 180
+  timeoutInMinutes: 240
 
   steps:
   # set up caches:

.github/workflows/cuda.yml

@@ -115,7 +115,7 @@ jobs:
         which nvcc || echo "nvcc not in PATH!"
 
         git clone https://github.com/AMReX-Codes/amrex.git ../amrex
-        cd ../amrex && git checkout --detach 68244ec91d118b5d4cc21f93376eaae8b56c51eb && cd -
+        cd ../amrex && git checkout --detach 2230caa24c7d4bd07edb08b54e9368f9c73eae6e && cd -
         make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_PSATD=TRUE USE_CCACHE=TRUE -j 4
 
         ccache -s

Docs/source/developers/particles.rst

@@ -36,17 +36,14 @@ A typical loop over particles reads:
       }
   }
 
-The innermost step ``[MY INNER LOOP]`` typically calls ``amrex::ParallelFor`` to perform operations on all particles in a portable way. For this reasons, the particle data needs to be converted in plain-old-data structures. The innermost loop in the code snippet above could look like:
+The innermost step ``[MY INNER LOOP]`` typically calls ``amrex::ParallelFor`` to perform operations on all particles in a portable way. The innermost loop in the code snippet above could look like:
 
 .. code-block:: cpp
 
-  // Get Array-Of-Struct particle data, also called data
-  // (x, y, z, id, cpu)
-  const auto& particles = pti.GetArrayOfStructs();
   // Get Struct-Of-Array particle data, also called attribs
-  // (ux, uy, uz, w, Exp, Ey, Ez, Bx, By, Bz)
+  // (x, y, z, ux, uy, uz, w)
   auto& attribs = pti.GetAttribs();
-  auto& Exp = attribs[PIdx::Ex];
+  auto& x = attribs[PIdx::x];
   // [...]
   // Number of particles in this box
   const long np = pti.numParticles();
@@ -66,7 +63,6 @@ On a loop over boxes in a ``MultiFab`` (``MFIter``), it can be useful to access
   const int tile_id = mfi.LocalTileIndex();
   // Get GPU-friendly arrays of particle data
   auto& ptile = GetParticles(lev)[std::make_pair(grid_id,tile_id)];
-  ParticleType* pp = particle_tile.GetArrayOfStructs()().data();
   // Only need attribs (i.e., SoA data)
   auto& soa = ptile.GetStructOfArrays();
   // As an example, let's get the ux momentum

Docs/source/usage/parameters.rst

@@ -2032,8 +2032,8 @@ Particle push, charge and current deposition, field gathering
 
      If ``algo.particle_pusher`` is not specified, ``boris`` is the default.
 
-* ``algo.particle_shape`` (`integer`; `1`, `2`, or `3`)
-    The order of the shape factors (splines) for the macro-particles along all spatial directions: `1` for linear, `2` for quadratic, `3` for cubic.
+* ``algo.particle_shape`` (`integer`; `1`, `2`, `3`, or `4`)
+    The order of the shape factors (splines) for the macro-particles along all spatial directions: `1` for linear, `2` for quadratic, `3` for cubic, `4` for quartic.
     Low-order shape factors result in faster simulations, but may lead to more noisy results.
     High-order shape factors are computationally more expensive, but may increase the overall accuracy of the results. For production runs it is generally safer to use high-order shape factors, such as cubic order.
 
@@ -2752,10 +2752,11 @@ 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 ``stepScraped``, which
-indicates the PIC iteration at which each particle was absorbed at the boundary,
-and a real attribute ``deltaTimeScraped``, which indicates the fraction of time between the time associated to the last step
-and the exact time when each particle hits the boundary. ``deltaTimeScraped`` is a real between ``0`` and ``dt``.
+In addition to their usual attributes, the saved particles have
+   an integer attribute ``stepScraped``, which indicates the PIC iteration at which each particle was absorbed at the boundary,
+   a real attribute ``deltaTimeScraped``, which indicates the time between the time associated to `stepScraped`
+   and the exact time when each particle hits the boundary.
+   3 real attributes ``nx``, ``ny``, ``nz``, which represents the three components of the normal to the boundary on the point of contact of the particles (not saved if they reach non-EB boundaries)
 
 ``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/langmuir/analysis_2d.py

@@ -38,6 +38,9 @@
 # Parse test name and check if Vay current deposition (algo.current_deposition=vay) is used
 vay_deposition = True if re.search( 'Vay_deposition', fn ) else False
 
+# Parse test name and check if particle_shape = 4 is used
+particle_shape_4 = True if re.search('particle_shape_4', fn) else False
+
 # Parameters (these parameters must match the parameters in `inputs.multi.rt`)
 epsilon = 0.01
 n = 4.e24
@@ -114,7 +117,11 @@ def get_theoretical_field( field, t ):
 fig.tight_layout()
 fig.savefig('Langmuir_multi_2d_analysis.png', dpi = 200)
 
-tolerance_rel = 0.05
+if particle_shape_4:
+# lower fidelity, due to smoothing
+    tolerance_rel = 0.07
+else:
+    tolerance_rel = 0.05
 
 print("error_rel    : " + str(error_rel))
 print("tolerance_rel: " + str(tolerance_rel))

Examples/Tests/point_of_contact_EB/analysis.py

@@ -26,35 +26,45 @@
 ts_scraping = OpenPMDTimeSeries('./diags/diag2/particles_at_eb/')
 
 it=ts_scraping.iterations
-step_scraped, delta, x, y, z=ts_scraping.get_particle( ['stepScraped','deltaTimeScraped','x','y','z'], species='electron', iteration=it )
-delta_reduced=delta[0]*1e10
+step_scraped, time_scraped, x, y, z, nx, ny, nz=ts_scraping.get_particle( ['stepScraped','timeScraped','x','y','z', 'nx', 'ny', 'nz'], 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
+nx_analytic=-0.99
+ny_analytic=0.13
+nz_analytic=0.0
 
 #result obtained by analysis of simulations
 step_ref=3
 delta_reduced_ref=0.59
 
 print('NUMERICAL coordinates of the point of contact:')
-print('step_scraped=%d, delta=%5.4f e-10, x=%5.4f, y=%5.4f, z=%5.4f' % (step_scraped[0],delta_reduced,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, nx=%5.4f, ny=%5.4f, nz=%5.4f' % (step_scraped[0],time_reduced,x[0], y[0], z[0], nx[0], ny[0], nz[0]))
 print('\n')
 print('ANALYTICAL coordinates of the point of contact:')
-print('step_scraped=%d, delta=%5.4f e-10, x=%5.4f, y=%5.4f, z=%5.4f' % (step_ref, delta_reduced_ref, 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, nx=%5.4f, ny=%5.4f, nz=%5.4f' % (step, time_reduced, x_analytic, y_analytic, z_analytic, nx_analytic, ny_analytic, nz_analytic))
 
 tolerance=0.001
-tolerance_t=0.01
+tolerance_t=0.003
+tolerance_n=0.01
 print("tolerance = "+ str(tolerance *100) + '%')
 print("tolerance for the time = "+ str(tolerance_t *100) + '%')
+print("tolerance for the normal components = "+ str(tolerance_n *100) + '%')
 
 diff_step=np.abs((step_scraped[0]-step_ref)/step_ref)
 diff_delta=np.abs((delta_reduced-delta_reduced_ref)/delta_reduced_ref)
 diff_x=np.abs((x[0]-x_analytic)/x_analytic)
 diff_y=np.abs((y[0]-y_analytic)/y_analytic)
+diff_nx=np.abs((nx[0]-nx_analytic)/nx_analytic)
+diff_ny=np.abs((ny[0]-ny_analytic)/ny_analytic)
 
 print("percentage error for x = %5.4f %%" %(diff_x *100))
 print("percentage error for y = %5.4f %%" %(diff_y *100))
+print("percentage error for nx = %5.2f %%" %(diff_nx *100))
+print("percentage error for ny = %5.2f %%" %(diff_ny *100))
+print("nz = %5.2f " %(nz[0]))
 
-assert (diff_x < tolerance) and (diff_y < tolerance) and (np.abs(z[0]) < 1e-8) and (diff_step < 1e-8) and (diff_delta < tolerance_t) , '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) and (diff_nx < tolerance_n) and (diff_ny < tolerance_n) and (np.abs(nz) < 1e-8) , 'Test point_of_contact did not pass'

Python/pywarpx/particle_containers.py

@@ -774,7 +774,9 @@ def get_particle_boundary_buffer(self, species_name, boundary, comp_name, level)
             comp_name      : str
                 The component of the array data that will be returned.
                 "x", "y", "z", "ux", "uy", "uz", "w"
-                "stepScraped","deltaTimeScraped", "nx", "ny", "nz"
+
+                "stepScraped","timeScraped",
+                if boundary='eb': "nx", "ny", "nz"
 
             level          : int
                 Which AMR level to retrieve scraped particle data from.

Regression/Checksum/benchmarks_json/Langmuir_multi_2d_psatd_Vay_deposition_particle_shape_4.json

@@ -0,0 +1,29 @@
+{
+  "electrons": {
+    "particle_momentum_x": 5.696397887862144e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 5.696397887862156e-20,
+    "particle_position_x": 0.6553600000001614,
+    "particle_position_y": 0.6553600000001614,
+    "particle_weight": 3200000000000000.5
+  },
+  "lev=0": {
+    "Ex": 3616987165668.129,
+    "Ey": 0.0,
+    "Ez": 3616987165667.9756,
+    "divE": 2.269072850514105e+18,
+    "jx": 1.0121499183289864e+16,
+    "jy": 0.0,
+    "jz": 1.0121499183289892e+16,
+    "part_per_cell": 131072.0,
+    "rho": 20090797.21028113
+  },
+  "positrons": {
+    "particle_momentum_x": 5.696397887862144e-20,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 5.696397887862156e-20,
+    "particle_position_x": 0.6553600000001614,
+    "particle_position_y": 0.6553600000001614,
+    "particle_weight": 3200000000000000.5
+  }
+}

Regression/WarpX-GPU-tests.ini

@@ -60,7 +60,7 @@ emailBody = Check https://ccse.lbl.gov/pub/GpuRegressionTesting/WarpX/ for more
 
 [AMReX]
 dir = /home/regtester/git/amrex/
-branch = 68244ec91d118b5d4cc21f93376eaae8b56c51eb
+branch = 2230caa24c7d4bd07edb08b54e9368f9c73eae6e
 
 [source]
 dir = /home/regtester/git/WarpX

Regression/WarpX-tests.ini

@@ -59,7 +59,7 @@ emailBody = Check https://ccse.lbl.gov/pub/RegressionTesting/WarpX/ for more det
 
 [AMReX]
 dir = /home/regtester/AMReX_RegTesting/amrex/
-branch = 68244ec91d118b5d4cc21f93376eaae8b56c51eb
+branch = 2230caa24c7d4bd07edb08b54e9368f9c73eae6e
 
 [source]
 dir = /home/regtester/AMReX_RegTesting/warpx
@@ -1434,6 +1434,25 @@ particleTypes = electrons positrons
 analysisRoutine = Examples/Tests/langmuir/analysis_2d.py
 analysisOutputImage = langmuir_multi_2d_analysis.png
 
+[Langmuir_multi_2d_psatd_Vay_deposition_particle_shape_4]
+buildDir = .
+inputFile = Examples/Tests/langmuir/inputs_2d
+runtime_params = algo.maxwell_solver=psatd amr.max_grid_size=128 algo.current_deposition=vay diag1.electrons.variables=w ux uy uz diag1.positrons.variables=w ux uy uz diag1.fields_to_plot = Ex Ey Ez jx jy jz part_per_cell rho divE warpx.cfl = 0.7071067811865475 algo.particle_shape=4
+dim = 2
+addToCompileString = USE_PSATD=TRUE
+cmakeSetupOpts = -DWarpX_DIMS=2 -DWarpX_PSATD=ON
+restartTest = 0
+useMPI = 1
+numprocs = 1
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+compareParticles = 1
+particleTypes = electrons positrons
+analysisRoutine = Examples/Tests/langmuir/analysis_2d.py
+analysisOutputImage = langmuir_multi_2d_analysis.png
+
 [Langmuir_multi_2d_psatd_Vay_deposition_nodal]
 buildDir = .
 inputFile = Examples/Tests/langmuir/inputs_2d

Source/EmbeddedBoundary/DistanceToEB.H

@@ -34,56 +34,95 @@ void normalize (amrex::RealVect& a) noexcept
                  a[2] *= inv_norm);
 }
 
+
+
+// This function calculates the normal vector using the nodal and cell-centered data.
+// i,j,k are the index of the nearest node to the left of the point at which we interpolate.
+// W are the interpolation weight for the left and right nodes (for the 0th component and 1st component respectively)
+// ic,jc,kc are the index of the nearest cell-center to the left of the point at which we interpolate.
 AMREX_GPU_HOST_DEVICE AMREX_INLINE
 amrex::RealVect interp_normal (int i, int j, int k, const amrex::Real W[AMREX_SPACEDIM][2],
+                               int ic, int jc, int kc, const amrex::Real Wc[AMREX_SPACEDIM][2],
                                amrex::Array4<const amrex::Real> const& phi,
                                amrex::GpuArray<amrex::Real,AMREX_SPACEDIM> const& dxi) noexcept
 {
+
 #if (defined WARPX_DIM_3D)
     amrex::RealVect normal{0.0, 0.0, 0.0};
+    for (int iic = 0; iic < 2; ++iic) {
+        for (int kk = 0; kk < 2; ++kk) {
+            for (int jj=0; jj< 2; ++jj) {
+                for (int ii = 0; ii < 2; ++ii) {
+                    int icstart = ic + iic;
+                    amrex::Real sign = (ii%2)*2. - 1.;
+                    int wccomp = static_cast<int>(iic%2);
+                    int w1comp = static_cast<int>(jj%2);
+                    int w2comp = static_cast<int>(kk%2);
+                    normal[0] += sign * phi(icstart + ii, j + jj, k + kk) * dxi[0] * Wc[0][wccomp] * W[1][w1comp] * W[2][w2comp];
+                }
+            }
+        }
+    }
+    for (int iic = 0; iic < 2; ++iic) {
+        for (int kk = 0; kk < 2; ++kk) {
+            for (int ii=0; ii< 2; ++ii) {
+                for (int jj = 0; jj < 2; ++jj) {
+                    int jcstart = jc + iic;
+                    amrex::Real sign = (jj%2)*2. - 1.;
+                    int wccomp = static_cast<int>(iic%2);
+                    int w1comp = static_cast<int>(ii%2);
+                    int w2comp = static_cast<int>(kk%2);
+                    normal[1] += sign * phi(i + ii, jcstart + jj, k + kk) * dxi[1] * W[0][w1comp] * Wc[1][wccomp] * W[2][w2comp];
+                }
+            }
+        }
+    }
+    for (int iic = 0; iic < 2; ++iic) {
+        for (int jj = 0; jj < 2; ++jj) {
+            for (int ii=0; ii< 2; ++ii) {
+                for (int kk = 0; kk < 2; ++kk) {
+                    int kcstart = kc + iic;
+                    amrex::Real sign = (kk%2)*2. - 1.;
+                    int wccomp = static_cast<int>(iic%2);
+                    int w1comp = static_cast<int>(ii%2);
+                    int w2comp = static_cast<int>(jj%2);
+                    normal[2] += sign * phi(i + ii, j + jj, kcstart + kk) * dxi[2] * W[0][w1comp] * W[1][w2comp] * Wc[2][wccomp];
+                }
+            }
+        }
+    }
 
-    normal[0] -= phi(i,   j  , k  ) * dxi[0] * W[1][0] * W[2][0];
-    normal[0] += phi(i+1, j  , k  ) * dxi[0] * W[1][0] * W[2][0];
-    normal[0] -= phi(i,   j+1, k  ) * dxi[0] * W[1][1] * W[2][0];
-    normal[0] += phi(i+1, j+1, k  ) * dxi[0] * W[1][1] * W[2][0];
-    normal[0] -= phi(i,   j  , k+1) * dxi[0] * W[1][0] * W[2][1];
-    normal[0] += phi(i+1, j  , k+1) * dxi[0] * W[1][0] * W[2][1];
-    normal[0] -= phi(i  , j+1, k+1) * dxi[0] * W[1][1] * W[2][1];
-    normal[0] += phi(i+1, j+1, k+1) * dxi[0] * W[1][1] * W[2][1];
-
-    normal[1] -= phi(i,   j  , k  ) * dxi[1] * W[0][0] * W[2][0];
-    normal[1] += phi(i  , j+1, k  ) * dxi[1] * W[0][0] * W[2][0];
-    normal[1] -= phi(i+1, j  , k  ) * dxi[1] * W[0][1] * W[2][0];
-    normal[1] += phi(i+1, j+1, k  ) * dxi[1] * W[0][1] * W[2][0];
-    normal[1] -= phi(i,   j  , k+1) * dxi[1] * W[0][0] * W[2][1];
-    normal[1] += phi(i  , j+1, k+1) * dxi[1] * W[0][0] * W[2][1];
-    normal[1] -= phi(i+1, j  , k+1) * dxi[1] * W[0][1] * W[2][1];
-    normal[1] += phi(i+1, j+1, k+1) * dxi[1] * W[0][1] * W[2][1];
-
-    normal[2] -= phi(i  , j  , k  ) * dxi[2] * W[0][0] * W[1][0];
-    normal[2] += phi(i  , j  , k+1) * dxi[2] * W[0][0] * W[1][0];
-    normal[2] -= phi(i+1, j  , k  ) * dxi[2] * W[0][1] * W[1][0];
-    normal[2] += phi(i+1, j  , k+1) * dxi[2] * W[0][1] * W[1][0];
-    normal[2] -= phi(i,   j+1, k  ) * dxi[2] * W[0][0] * W[1][1];
-    normal[2] += phi(i  , j+1, k+1) * dxi[2] * W[0][0] * W[1][1];
-    normal[2] -= phi(i+1, j+1, k  ) * dxi[2] * W[0][1] * W[1][1];
-    normal[2] += phi(i+1, j+1, k+1) * dxi[2] * W[0][1] * W[1][1];
 #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
     amrex::RealVect normal{0.0, 0.0};
+    for (int iic = 0; iic < 2; ++iic) {
+        for (int jj=0; jj< 2; ++jj) {
+            for (int ii = 0; ii < 2; ++ii) {
+                int icstart = ic + iic;
+                amrex::Real sign = (ii%2)*2. - 1.;
+                int wccomp = static_cast<int>(iic%2);
+                int w1comp = static_cast<int>(jj%2);
+                normal[0] += sign * phi(icstart + ii, j + jj, k) * dxi[0] * Wc[0][wccomp] * W[1][w1comp];
+            }
+        }
+    }
+    for (int iic = 0; iic < 2; ++iic) {
+        for (int ii=0; ii< 2; ++ii) {
+            for (int jj = 0; jj < 2; ++jj) {
+                int jcstart = jc + iic;
+                amrex::Real sign = (jj%2)*2. - 1.;
+                int wccomp = static_cast<int>(iic%2);
+                int w1comp = static_cast<int>(ii%2);
+                normal[1] += sign * phi(i + ii, jcstart + jj, k) * dxi[1] * W[0][w1comp] * Wc[1][wccomp];
+            }
+        }
+    }
+    amrex::ignore_unused(kc);
 
-    normal[0] -= phi(i,   j  , k) * dxi[0] * W[1][0];
-    normal[0] += phi(i+1, j  , k) * dxi[0] * W[1][0];
-    normal[0] -= phi(i,   j+1, k) * dxi[0] * W[1][1];
-    normal[0] += phi(i+1, j+1, k) * dxi[0] * W[1][1];
-
-    normal[1] -= phi(i,   j  , k) * dxi[1] * W[0][0];
-    normal[1] += phi(i  , j+1, k) * dxi[1] * W[0][0];
-    normal[1] -= phi(i+1, j  , k) * dxi[1] * W[0][1];
-    normal[1] += phi(i+1, j+1, k) * dxi[1] * W[0][1];
 #else
+    amrex::ignore_unused(i, j, k, ic, jc, kc, W, Wc, phi, dxi);
     amrex::RealVect normal{0.0, 0.0};
-    amrex::ignore_unused(i, j, k, W, phi, dxi);
     WARPX_ABORT_WITH_MESSAGE("Error: interp_distance not yet implemented in 1D");
+
 #endif
     return normal;
 }

Source/EmbeddedBoundary/ParticleScraper.H

@@ -10,6 +10,7 @@
 #include "EmbeddedBoundary/DistanceToEB.H"
 #include "Particles/Pusher/GetAndSetPosition.H"
 
+
 #include <ablastr/particles/NodalFieldGather.H>
 
 #include <AMReX.H>
@@ -182,16 +183,17 @@ scrapeParticles (PC& pc, const amrex::Vector<const amrex::MultiFab*>& distance_t
 
                 int i, j, k;
                 amrex::Real W[AMREX_SPACEDIM][2];
-                ablastr::particles::compute_weights_nodal(xp, yp, zp, plo, dxi, i, j, k, W);
-
+                ablastr::particles::compute_weights(xp, yp, zp, plo, dxi, i, j, k, W);
                 amrex::Real phi_value  = ablastr::particles::interp_field_nodal(i, j, k, W, phi);
 
                 if (phi_value < 0.0)
                 {
-
-                    amrex::RealVect normal = DistanceToEB::interp_normal(i, j, k, W, phi, dxi);
+                    int ic, jc, kc; // Cell-centered indices
+                    int nodal;
+                    amrex::Real Wc[AMREX_SPACEDIM][2]; // Cell-centered weights
+                    ablastr::particles::compute_weights(xp, yp, zp, plo, dxi, ic, jc, kc, Wc, nodal=0);
+                    amrex::RealVect normal = DistanceToEB::interp_normal(i, j, k, W, ic, jc, kc, Wc, phi, dxi);
                     DistanceToEB::normalize(normal);
-
                     amrex::RealVect pos;
 #if (defined WARPX_DIM_3D)
                     pos[0] = xp;