Merge branch 'develop'

benchmarks/tstAM_04d_laser_propagation_Order1.py

@@ -0,0 +1,98 @@
+# ----------------------------------------------------------------------------------------
+# 					SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI
+# ----------------------------------------------------------------------------------------
+
+import math
+dx = 0.251327
+dtrans = 1.96349
+dt = 0.96 * dx
+nx =  960
+ntrans = 256
+Lx = nx * dx
+Ltrans = ntrans * dtrans
+npatch_x = 64
+npatch_trans =32
+Nit = 2000
+
+
+Main(
+    geometry = "AMcylindrical",
+    number_of_AM=2,
+    interpolation_order = 1,
+    timestep = dt,
+    simulation_time = dt*Nit,
+    cell_length  = [dx, dtrans],
+    grid_length = [ Lx,  Ltrans],
+    number_of_patches = [npatch_x, npatch_trans],
+    cluster_width = 5,
+    EM_boundary_conditions = [
+        ["silver-muller","silver-muller"],
+        ["PML","PML"],
+    ],
+    number_of_pml_cells = [[0,0], [10,10]],
+    solve_poisson = False,
+    print_every = 100
+)
+
+MovingWindow(
+    time_start = Main.grid_length[0] - 50*dx, #Leaves 2 patches untouched, in front of the laser.
+    velocity_x = 0.996995486
+)
+
+ne = 0.0045
+begin_upramp = 10.  #Density is 0 before that and up ramp starts.
+Lupramp = 100. #Length of the upramp 
+Lplateau = 15707.  #Length of the plateau 
+Ldownramp = 2356.19 #Length of the down ramp
+xplateau = begin_upramp + Lupramp # Start of the plateau
+begin_downramp = xplateau + Lplateau # Beginning of the output ramp. 
+finish = begin_downramp + Ldownramp # End of plasma
+
+g = polygonal(xpoints=[begin_upramp, xplateau, begin_downramp, finish], xvalues=[0, ne, ne, 0.])
+
+def my_profile(x,y):
+    return g(x,y)
+
+Species( 
+    name = "electron",
+    position_initialization = "regular",
+    momentum_initialization = "cold",
+    ionization_model = "none",
+    particles_per_cell = 30,
+    c_part_max = 1.0,
+    mass = 1.0,
+    charge = -1.0,
+    charge_density = my_profile,  # Here absolute value of the charge is 1 so charge_density = nb_density
+    mean_velocity = [0., 0., 0.],
+    time_frozen = 0.0,
+    boundary_conditions = [
+    	["remove", "remove"],
+    	["reflective", "remove"],
+    ],
+)
+
+laser_fwhm = 82. 
+LaserGaussianAM(
+    box_side         = "xmin",
+    a0              = 2.,
+    focus           = [10.],  
+    waist           = 120.,
+    time_envelope   = tgaussian(center=2**0.5*laser_fwhm, fwhm=laser_fwhm)
+)
+
+DiagProbe(
+	every = 1000,
+	origin = [0., 2*dtrans, 0.],
+	corners = [
+              [Main.grid_length[0], 2*dtrans, 0.]
+                  ],
+	number = [nx],
+)
+
+DiagPerformances(
+	every = 1000,
+)
+
+DiagScalar(
+	every = 20
+)

doc/Sphinx/Overview/releases.rst

@@ -16,13 +16,27 @@ Get Smilei
 
 You can find older, `unsupported versions here <https://github.com/SmileiPIC/Smilei/releases>`_
 
-..
-.. ----
+----
+
+.. _latestVersion:
+
+Changes made in the repository (not released)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+* **Features**:
+
+  * Prescribed fields in AM geometry.
+  * Particle reflective boundary conditions at Rmax in AM geometry.
+  * 1st order Ruyten shape function in AM geometry.
+
+* **Bug fixes**:
 
-.. .. _latestVersion:
+  * Tunnel ionization was wrong in some cases for high atomic numbers.
+  * Custom functions in ``ParticleBinning`` crashed with python 3.12.
+  * Species-specific diagnostics in AM geometry with vectorization.
+  * Happi's ``average`` argument would sometimes be missing the last bin.
+  * 1D projector on GPU without diagnostics
 
-.. Changes made in the repository (not released)
-.. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 ----
 

doc/Sphinx/Understand/GPU_offloading.rst

@@ -17,9 +17,11 @@ the announced exaflopic supercomputers will include GPUs.
 * Currently supported features:
 
   * Both AMD's GPUs and Nvidia's GPUs are supported
-  * Cartesian geometry in 2D and in 3D
+  * Cartesian geometry in 1D, 2D and in 3D , for order 2
   * Diagnostics: Field, Probes, Scalar, ParticleBinning, TrackParticles
-  * Moving Window.
+  * Moving Window
 
 * A few key features remain to be implemented (AM geometry, ionization, PML, envelope,
-  additional physics), but the fundamentals of the code are ported.
+  additional physics), but the fundamentals of the code are ported.
+
+* Short term roadmap We are currently working on porting on GPU the following features: AM geometry (order 2) and collisions

doc/Sphinx/Use/GPU_version.rst

@@ -0,0 +1,23 @@
+GPU version of SMILEI
+-------
+
+
+This page contains the links of this documentation to compile and run SMILEI on GPU clusters, as well as the list of the features that are supported or in development
+
+* :doc:`Introduction</Overview/highlights>` #Improved performance using GPU offloading
+
+* :doc:`GPU offloading , supported features, guidelines and roadmap</Understand/GPU_offloading>`
+
+* :doc:`Intallation and compilation</Use/install_linux_GPU>`
+
+* :doc:`Compilation for GPU-accelerated nodes</Use/installation>` #compilation-for-gpu-accelerated-nodes
+
+* :doc:`Running on GPU-equiped nodes</Use/run>` #running-on-gpu-equiped-nodes
+
+----
+
+Known issues
+^^^^^^^^^^^^
+
+2D and 3D runs may crash with A2000 & A6000 GPUs (used in laptops and worstations respectively, 
+they are not 'production GPUs' which are designed for 64 bits floating point operations )

doc/Sphinx/Use/installation.rst

@@ -185,8 +185,6 @@ Typically ``CXXFLAGS += -ta=tesla:cc80`` for ``nvhpc`` <23.4 and
 
 .. warning::
 
-  * We are aware of issues with CUDA >12.0, fixes are being tested but are not deployed yet.
-    We recommend CUDA 11.x at the moment.
   * The hdf5 module should be compiled with the nvidia/cray compiler;
     openmpi as well, but depending on the nvhpc module it might not be needed as it can be included in the nvhpc module.
 

doc/Sphinx/Use/namelist.rst

@@ -143,9 +143,13 @@ The block ``Main`` is **mandatory** and has the following syntax::
 
   Interpolation order, defines particle shape function:
 
+  * ``1``  : 2 points stencil in r with Ruyten correction, 3 points stencil in x. Supported only in AM geometry.
   * ``2``  : 3 points stencil, supported in all configurations.
   * ``4``  : 5 points stencil, not supported in vectorized 2D geometry.
 
+  The Ruyten correction is the scheme described bu equation 4.2 in `this paper <https://www.sciencedirect.com/science/article/abs/pii/S0021999183710703>`_ .
+  It allows for a more accurate description on axis at the cost of a higher statistic noise so it often requires the use of more macro-particles.
+
 .. py:data:: interpolator
 
   :default: ``"momentum-conserving"``
@@ -2009,8 +2013,8 @@ at the beginning of the simulation using the ``ExternalField`` block::
 
 .. py:data:: field
 
-  Field name in Cartesian geometries: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx"``, ``"By"``, ``"Bz"``, ``"Bx_m"``, ``"By_m"``, ``"Bz_m"``
-  Field name in AM geometry: ``"El"``, ``"Er"``, ``"Et"``, ``"Bl"``, ``"Br"``, ``"Bt"``, ``"Bl_m"``, ``"Br_m"``, ``"Bt_m"``, ``"A"``, ``"A0"`` .
+  Field names in Cartesian geometries: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx"``, ``"By"``, ``"Bz"``, ``"Bx_m"``, ``"By_m"``, ``"Bz_m"``.
+  Field names in AM geometry: ``"El_mode_m"``, ``"Er_mode_m"``, ``"Et_mode_m"``, ``"Bl_mode_m"``, ``"Br_mode_m"``, ``"Bt_mode_m"``, ``"Bl_m_mode_m"``, ``"Br_m_mode_m"``, ``"Bt_m_mode_m"``, ``"A_mode_1"``, ``"A0_mode_1"`` .
 
 .. py:data:: profile
 
@@ -2055,13 +2059,25 @@ This feature is accessible using the ``PrescribedField`` block::
 
 .. py:data:: field
 
-  Field name: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``.
+  Field names in Cartesian geometries: ``"Ex"``, ``"Ey"``, ``"Ez"``, ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``.
+  Field names in AM geometry: ``"El_mode_m"``, ``"Er_mode_m"``, ``"Et_mode_m"``, ``"Bl_m_mode_m"``, ``"Br_m_mode_m"`` or ``"Bt_m_mode_m"``.
 
 .. warning::
 
   When prescribing a magnetic field, always use the time-centered fields ``"Bx_m"``, ``"By_m"`` or ``"Bz_m"``.
   These fields are those used in the particle pusher, and are defined at integer time-steps.
 
+.. warning::
+
+  When prescribing a field in AM geometry, the mode "m" must be specified explicitly in the name of the field and the profile
+  must return a complex value.
+
+.. warning::
+
+  ``PrescribedFields`` are not visible in the ``Field`` diagnostic, 
+  but can be visualised through ``Probes`` and with the fields attributes of ``TrackParticles`` 
+  (since they sample the total field acting on the macro-particles).
+
 .. py:data:: profile
 
   :type: float or :doc:`profile <profiles>`
@@ -3641,4 +3657,4 @@ namelist. They should not be re-defined by the user!
 .. note::
 
   These variables can be access during ``happi`` post-processing, e.g.
-  ``S.namelist.smilei_mpi_size``.
+  ``S.namelist.smilei_mpi_size``.

doc/Sphinx/use.rst

@@ -12,3 +12,4 @@ Use
    Use/post-processing
    Use/contribute
    Use/troubleshoot
+   Use/GPU_version

happi/_Diagnostics/Diagnostic.py

@@ -578,8 +578,8 @@ def _selectRange(self, portion, meshpoints, axisname, axisunits, operation, edge
 					info = operation+" for "+axisname+" from "+axismin+" to "+axismax
 				else:
 					info = operation+" for "+axisname+" from "+str(meshpoints[indices[0]])+" to "+str(meshpoints[indices[-1]])+" "+axisunits
-				selection = slice(indices[0],indices[-1])
-				finalShape = indices[-1] - indices[0]
+				selection = slice(indices[0], indices[-1] + 1)
+				finalShape = indices[-1] - indices[0] + 1
 		return info, selection, finalShape
 
 	# Method to prepare some data before plotting

src/Diagnostic/DiagnosticNewParticles.cpp

@@ -93,8 +93,8 @@ void DiagnosticNewParticles::openFile( Params &params, SmileiMPI *smpi )
         }
     }
     if( write_weight_ ) {
-        loc_weight_ = newDataset( species_group, "weight", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_DENSITY );
-        openPMD_->writeRecordAttributes( *loc_weight_, SMILEI_UNIT_DENSITY );
+        loc_weight_ = newDataset( species_group, "weight", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_WEIGHT );
+        openPMD_->writeRecordAttributes( *loc_weight_, SMILEI_UNIT_WEIGHT );
     }
     if( write_chi_ ) {
         loc_chi_ = newDataset( species_group, "chi", H5T_NATIVE_DOUBLE, file_space, SMILEI_UNIT_NONE );

src/Diagnostic/DiagnosticParticleList.cpp

@@ -299,7 +299,7 @@ void DiagnosticParticleList::run( SmileiMPI *smpi, VectorPatch &vecPatches, int
     if( write_weight_ ) {
         fill_buffer( vecPatches, iprop, data_double );
         #pragma omp master
-        write_scalar_double( loc_weight_, "weight", data_double[0], file_space, mem_space, SMILEI_UNIT_DENSITY );
+        write_scalar_double( loc_weight_, "weight", data_double[0], file_space, mem_space, SMILEI_UNIT_WEIGHT );
     }
 
     iprop++;

src/Diagnostic/Histogram.h

@@ -578,25 +578,24 @@ class HistogramAxis_user_function : public HistogramAxis
 private:
     void calculate_locations( Species *s, double *array, int *index, unsigned int npart, SimWindow * )
     {
-        #pragma omp critical
+        SMILEI_PY_ACQUIRE_GIL
+        // Expose particle data as numpy arrays
+        particleData.resize( npart );
+        particleData.set( s->particles );
+        // run the function
+        PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL );
+        particleData.clear();
+        // Copy the result to "array"
+        double *arr = ( double * ) PyArray_GETPTR1( ret, 0 );
+        for( unsigned int ipart = 0 ; ipart < npart ; ipart++ )
         {
-            // Expose particle data as numpy arrays
-            particleData.resize( npart );
-            particleData.set( s->particles );
-            // run the function
-            PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL );
-            particleData.clear();
-            // Copy the result to "array"
-            double *arr = ( double * ) PyArray_GETPTR1( ret, 0 );
-            for( unsigned int ipart = 0 ; ipart < npart ; ipart++ )
-            {
-                if( index[ipart]<0 ) {
-                    continue;
-                }
-                array[ipart] = arr[ipart];
+            if( index[ipart]<0 ) {
+                continue;
             }
-            Py_DECREF( ret );
+            array[ipart] = arr[ipart];
         }
+        Py_DECREF( ret );
+        SMILEI_PY_RELEASE_GIL
     };
 
     PyObject *function;
@@ -1167,26 +1166,25 @@ class Histogram_user_function : public Histogram
     void valuate( Species *s, double *array, int *index )
     {
         unsigned int npart = s->getNbrOfParticles();
-        #pragma omp critical
-        {
-            // Expose particle data as numpy arrays
-            particleData.resize( npart );
-            particleData.set( s->particles );
-            // run the function
-            PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL );
-            particleData.clear();
-            // Copy the result to "array"
-            double *arr = ( double * ) PyArray_GETPTR1( ret, 0 );
-            for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) {
-                if( index[ipart]<0 ) {
-                    continue;
-                }
-                array[ipart] = arr[ipart];
+        SMILEI_PY_ACQUIRE_GIL
+        // Expose particle data as numpy arrays
+        particleData.resize( npart );
+        particleData.set( s->particles );
+        // run the function
+        PyArrayObject *ret = ( PyArrayObject * )PyObject_CallFunctionObjArgs( function, particleData.get(), NULL );
+        particleData.clear();
+        // Copy the result to "array"
+        double *arr = ( double * ) PyArray_GETPTR1( ret, 0 );
+        for( unsigned int ipart = 0 ; ipart < npart ; ipart++ ) {
+            if( index[ipart]<0 ) {
+                continue;
             }
-            Py_DECREF( ret );
+            array[ipart] = arr[ipart];
         }
+        Py_DECREF( ret );
+        SMILEI_PY_RELEASE_GIL
     };
-
+    
     PyObject *function;
     ParticleData particleData;
 };