Merge branch 'develop' into develop

doc/Sphinx/Overview/material.rst

@@ -71,8 +71,6 @@ As of November 2021, 90 papers have been published covering a broad range of top
      `Modeling of a Liquid Leaf Target TNSA Experiment Using Particle-In-Cell Simulations and Deep Learning`,
      `Laser and Particle Beams, 2868112 (2023) <https://doi.org/10.1155/2023/2868112>`_
      
-
-       
 .. [Paschke_Bruehl2023]
 
     F. Paschke-Bruehl, M. Banjafar, M. Garten, L. G. Huang, B. E. Marré, M. Nakatsutsumi, L. Randolph, T. E. Cowan, U. Schramm and T. Kluge,
@@ -371,8 +369,8 @@ As of November 2021, 90 papers have been published covering a broad range of top
   `In International Conference on High Performance Computing in Asia-Pacific Region Workshops (HPCAsia 2022 Workshop). 
   Association for Computing Machinery, New York, NY, USA, 40–48. (2022) <http://doi.org/10.1145/3503470.3503475>`_
 
-  
-  .. [Tomassini2021]
+
+.. [Tomassini2021]
 
     Paolo Tomassini, Francesco Massimo, Luca Labate and Leonida A. Gizzi,
     `Accurate electron beam phase-space theory for ionization-injection schemes driven by laser pulses`,

doc/Sphinx/Overview/releases.rst

@@ -23,7 +23,11 @@ You can find older, `unsupported versions here <https://github.com/SmileiPIC/Smi
 Changes made in the repository (not released)
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-...
+* Happi:
+
+  * New arguments ``xoffset`` and ``yoffset`` to shift plot coordinates.
+  * Changed coordinate reference for 2D probe in 3D or AM geometry
+    (zero is the box origin projected orthogonally on the probe plane).
 
 ----
 

doc/Sphinx/Use/namelist.rst

@@ -160,8 +160,9 @@ The block ``Main`` is **mandatory** and has the following syntax::
              number_of_cells
 
   A list of numbers: size of the simulation box for each dimension of the simulation.
-   * Either ``grid_length``, the simulation length in each direction in units of :math:`L_r`,
-   * or ``number_of_cells``, the number of cells in each direction.
+
+  * Either ``grid_length``, the simulation length in each direction in units of :math:`L_r`,
+  * or ``number_of_cells``, the number of cells in each direction.
 
 
 .. py:data:: cell_length
@@ -173,20 +174,23 @@ The block ``Main`` is **mandatory** and has the following syntax::
              number_of_timesteps
 
   Duration of the simulation.
-    * Either ``simulation_time``, the simulation duration in units of :math:`T_r`,
-    * or ``number_of_timesteps``, the total number of timesteps.
+
+  * Either ``simulation_time``, the simulation duration in units of :math:`T_r`,
+  * or ``number_of_timesteps``, the total number of timesteps.
 
 
 .. py:data:: timestep
              timestep_over_CFL
 
   Duration of one timestep.
-    * Either ``timestep``, in units of :math:`T_r`,
-    * or ``timestep_over_CFL``, in units of the *Courant–Friedrichs–Lewy* (CFL) time.
+
+  * Either ``timestep``, in units of :math:`T_r`,
+  * or ``timestep_over_CFL``, in units of the *Courant–Friedrichs–Lewy* (CFL) time.
 
 .. py:data:: gpu_computing
 
    :default: ``False``
+
    Activates GPU acceleration if set to True
 
 .. py:data:: number_of_patches
@@ -430,11 +434,11 @@ The block ``Main`` is **mandatory** and has the following syntax::
   The number of azimuthal modes used for the relativistic field initialization in ``"AMcylindrical"`` geometry.
   Note that this number must be lower or equal to the number of modes of the simulation.
 
-  .. py:data:: use_BTIS3_interpolation
+.. py:data:: use_BTIS3_interpolation
 
-    :default: ``False``
+  :default: ``False``
 
-    If ``True``, the B-translated interpolation scheme 3 (or B-TIS3) described in :doc:`/Understand/algorithms` is used.
+  If ``True``, the B-translated interpolation scheme 3 (or B-TIS3) described in :doc:`/Understand/algorithms` is used.
 
 .. py:data:: custom_oversize
 
@@ -920,6 +924,7 @@ Each species has to be defined in a ``Species`` block::
 .. py:data:: temperature
 
   :type: a list of 3 floats or :doc:`profiles <profiles>`
+  :default: ``1e-10``
 
   The initial temperature of the particles, in units of :math:`m_ec^2`.
 

doc/Sphinx/Use/post-processing.rst

@@ -557,6 +557,9 @@ to manipulate the plotting options:
 
 * ``xmin``, ``xmax``, ``ymin``, ``ymax``: axes limits.
 * ``xfactor``, ``yfactor``: factors to rescale axes.
+* ``xoffset``, ``yoffset``: numerical values to offset the
+  coordinates. These values must be given in the original (normalized)
+  units, i.e. not acounting for the factors above or for unit conversion.
 * ``title``: a string that replaces the plot title (or the y-label in a 1D plot).
   The current simulation time can be included with the placeholders ``{time}`` and
   ``{time_units}``, together with formatting instructions conforming to

doc/Sphinx/Use/run.rst

@@ -128,7 +128,7 @@ administrators.
 ----
 
 Running on GPU-equiped nodes
-^^^^^^^^^^^^^^^^^^^^^^^^^
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 On a supercomputer equipped with GPUs it is necessary to use a binding script.
 Here are two examples: 

happi/_Diagnostics/Diagnostic.py

@@ -131,10 +131,12 @@ def limits(self):
 		"""
 		assert self.dim <= 2, "Method limits() may only be used in 1D or 2D"
 		self._prepare1()
-		l = []
 		factor = [self._xfactor, self._yfactor]
-		for i in range(self.dim):
-			l.append([self._centers[i].min()*factor[i], self._centers[i].max()*factor[i]])
+		offset = [self._xoffset, self._yoffset]
+		l = [[
+				(offset[i] + self._centers[i].min())*factor[i], 
+				(offset[i] + self._centers[i].max())*factor[i]
+			] for i in range(self.dim)]
 		return l
 
 	# Method to print info on this diag
@@ -213,12 +215,15 @@ def getAxis(self, axis, timestep=0):
 		except Exception as e: return []
 		if   axis_index == 0:
 			factor = (self.options.xfactor or 1.) * self.units.xcoeff
+			offset = self.options.xoffset or 0.
 		elif axis_index == 1:
 			factor = (self.options.yfactor or 1.) * self.units.ycoeff
+			offset = self.options.yoffset or 0.
 		else:
 			factor, _ = self.units._convert(self._units[axis_index], None)
+			offset = 0.
 		axis = self._getCenters(axis_index, timestep)
-		return factor * axis
+		return factor * ( offset + axis )
 
 	# Method to obtain the data and the axes
 	def get(self, timestep=None):
@@ -355,8 +360,8 @@ def streak(self, saveAs=None, axes=None, dpi=200, **kwargs):
 		A = self._np.double([self._dataAtTime(t) for t in self._timesteps])
 		# Plot
 		ax.cla()
-		xmin = self._xfactor*self._centers[0][0]
-		xmax = self._xfactor*self._centers[0][-1]
+		xmin = self._xfactor*( self._xoffset + self._centers[0][0] )
+		xmax = self._xfactor*( self._xoffset + self._centers[0][-1] )
 		extent = [xmin, xmax, self._yfactor*self._timesteps[0], self._yfactor*self._timesteps[-1]]
 		if self._log[0]: extent[0:2] = [self._np.log10(xmin), self._np.log10(xmax)]
 		im = ax.imshow(self._np.flipud(A), vmin = self.options.vmin, vmax = self.options.vmax, extent=extent, **self.options.image)
@@ -609,6 +614,9 @@ def _prepare1(self):
 		self._yfactor = (self.options.yfactor or 1.) * self.units.ycoeff
 		self._vfactor = (self.options.vfactor or 1.) * self.units.vcoeff
 		self._tfactor = (self.options.xfactor or 1.) * self.units.tcoeff * self.timestep
+		# prepare the offsets
+		self._xoffset = self.options.xoffset or 0.
+		self._yoffset = self.options.yoffset or 0.
 	def _prepare2(self):
 		# prepare the animating function
 		if not self._plotOnAxes:
@@ -640,19 +648,7 @@ def _prepare2(self):
 			if self.options.yfactor: self._ylabel += "/"+str(self.options.yfactor)
 			self._ylabel = self._label[1] + " (" + self._ylabel + ")"
 			if self._log[1]: self._ylabel = "Log[ "+self._ylabel+" ]"
-			# prepare extent for 2d plots
-			self._extent = [
-				self._xfactor*self._centers[0][0],
-				self._xfactor*self._centers[0][-1],
-				self._yfactor*self._centers[1][0],
-				self._yfactor*self._centers[1][-1]
-			]
-			if self._log[0]:
-				self._extent[0] = self._np.log10(self._extent[0])
-				self._extent[1] = self._np.log10(self._extent[1])
-			if self._log[1]:
-				self._extent[2] = self._np.log10(self._extent[2])
-				self._extent[3] = self._np.log10(self._extent[3])
+			self._prepareExtent()
 		# prepare title
 		self._vlabel = ""
 		if self.units.vname: self._vlabel += " (" + self.units.vname + ")"
@@ -694,6 +690,21 @@ def _prepare3(self):
 
 	def _prepare4(self): pass
 
+	def _prepareExtent(self):
+		# prepare extent for 2d plots
+		self._extent = [
+			self._xfactor*( self._xoffset + self._centers[0][0]  ),
+			self._xfactor*( self._xoffset + self._centers[0][-1] ),
+			self._yfactor*( self._yoffset + self._centers[1][0]  ),
+			self._yfactor*( self._yoffset + self._centers[1][-1] )
+		]
+		if self._log[0]:
+			self._extent[0] = self._np.log10(self._extent[0])
+			self._extent[1] = self._np.log10(self._extent[1])
+		if self._log[1]:
+			self._extent[2] = self._np.log10(self._extent[2])
+			self._extent[3] = self._np.log10(self._extent[3])
+
 	# Method to set limits to a plot
 	def _setLimits(self, ax, xmin=None, xmax=None, ymin=None, ymax=None):
 		ax.autoscale(tight=True)

happi/_Diagnostics/Field.py

@@ -342,7 +342,7 @@ def _getDataAtTime(self, t):
 
 		# Handle moving window
 		if self.moving and "x_moved" in h5item.attrs and 'x' in self._type:
-			self._xoffset = h5item.attrs["x_moved"]
+			self._xoffset = h5item.attrs["x_moved"] + (self.options.xoffset or 0.)
 			if self.dim>1 and hasattr(self,"_extent"):
 				self._extent[0] = self._xfactor*(self._xoffset + self._centers[0][ 0])
 				self._extent[1] = self._xfactor*(self._xoffset + self._centers[0][-1])
@@ -385,7 +385,7 @@ def _theta_getDataAtTime(self, t):
 
 		# Handle moving window
 		if self.moving and "x_moved" in h5item.attrs and 'x' in self._type:
-			self._xoffset = h5item.attrs["x_moved"]
+			self._xoffset = h5item.attrs["x_moved"] + (self.options.xoffset or 0.)
 			if self.dim>1 and hasattr(self,"_extent"):
 				self._extent[0] = self._xfactor*(self._xoffset + self._centers[0][ 0])
 				self._extent[1] = self._xfactor*(self._xoffset + self._centers[0][-1])
@@ -442,7 +442,7 @@ def _build3d_getDataAtTime(self, t):
 
 		# Handle moving window
 		if self.moving and "x_moved" in h5item.attrs and 'x' in self._type:
-			self._xoffset = h5item.attrs["x_moved"]
+			self._xoffset = h5item.attrs["x_moved"] + (self.options.xoffset or 0.)
 			if self.dim>1 and hasattr(self,"_extent"):
 				self._extent[0] = self._xfactor*(self._xoffset + self._centers[0][ 0])
 				self._extent[1] = self._xfactor*(self._xoffset + self._centers[0][-1])

happi/_Diagnostics/ParticleBinning.py

@@ -498,19 +498,7 @@ def _getDataAtTime(self, t):
 		if self.auto_axes:
 			self._updateAxes(t)
 			if len(self._shape) > 1:
-				# prepare extent for 2d plots
-				self._extent = [
-					self._xfactor*self._centers[0][0],
-					self._xfactor*self._centers[0][-1],
-					self._yfactor*self._centers[1][0],
-					self._yfactor*self._centers[1][-1]
-				]
-				if self._log[0]:
-					self._extent[0] = self._np.log10(self._extent[0])
-					self._extent[1] = self._np.log10(self._extent[1])
-				if self._log[1]:
-					self._extent[2] = self._np.log10(self._extent[2])
-					self._extent[3] = self._np.log10(self._extent[3])
+				self._prepareExtent()
 		# Get arrays from all requested diagnostics
 		A = {}
 		for d in self._diags:

happi/_Diagnostics/Performances.py

@@ -442,12 +442,11 @@ def toVTK(self,numberOfPieces=1,axis_quantity="patch"):
 
 			if self._verbose: print("Successfully exported to VTK, folder='"+self._exportDir)
 
-	def _prepare4(self):
-		if self._mode == "map" and self._ndim_fields==2:
-			self._extent = [
-				0., self._xfactor*self._number_of_patches[0]*self._patch_length[0],
-				0., self._yfactor*self._number_of_patches[1]*self._patch_length[1]
-			]
+	def _prepareExtent(self):
+		self._extent = [
+			self._xfactor*self._xoffset, self._xfactor*( self._xoffset + self._number_of_patches[0]*self._patch_length[0] ),
+			self._yfactor*self._yoffset, self._yfactor*( self._yoffset + self._number_of_patches[1]*self._patch_length[1] )
+		]
 
 	def _calculateMPIcontours_2D(self):
 		# Add lines to visualize MPI contours
@@ -462,9 +461,9 @@ def _calculateMPIcontours_2D(self):
 			vlines_i += [ self._np.full((len(j)//2), i, dtype=self._np.uint32) ]
 			vlines_jmin += [ j[ 0::2 ] ]
 			vlines_jmax += [ j[ 1::2 ] ]
-		vlines_i    = self._np.concatenate( vlines_i    )*self._xfactor*self._patch_length[0]
-		vlines_jmin = self._np.concatenate( vlines_jmin )*self._yfactor*self._patch_length[1]
-		vlines_jmax = self._np.concatenate( vlines_jmax )*self._yfactor*self._patch_length[1]
+		vlines_i    = (self._xoffset + self._np.concatenate( vlines_i    )*self._patch_length[0])*self._xfactor
+		vlines_jmin = (self._yoffset + self._np.concatenate( vlines_jmin )*self._patch_length[1])*self._yfactor
+		vlines_jmax = (self._yoffset + self._np.concatenate( vlines_jmax )*self._patch_length[1])*self._yfactor
 
 		# Horizontal lines
 		hlines_j    = []
@@ -477,9 +476,9 @@ def _calculateMPIcontours_2D(self):
 			hlines_j += [ self._np.full((len(i)//2), j, dtype=self._np.uint32) ]
 			hlines_imin += [ i[ 0::2 ] ]
 			hlines_imax += [ i[ 1::2 ] ]
-		hlines_j    = self._np.concatenate( hlines_j    )*self._yfactor*self._patch_length[1]
-		hlines_imin = self._np.concatenate( hlines_imin )*self._xfactor*self._patch_length[0]
-		hlines_imax = self._np.concatenate( hlines_imax )*self._xfactor*self._patch_length[0]
+		hlines_j    = (self._yoffset + self._np.concatenate( hlines_j    )*self._patch_length[1])*self._yfactor
+		hlines_imin = (self._xoffset + self._np.concatenate( hlines_imin )*self._patch_length[0])*self._xfactor
+		hlines_imax = (self._xoffset + self._np.concatenate( hlines_imax )*self._patch_length[0])*self._xfactor
 
 		return vlines_i, vlines_jmin, vlines_jmax, hlines_j, hlines_imin, hlines_imax
 

happi/_Diagnostics/Probe.py

@@ -162,29 +162,30 @@ def _init(self, requestedProbe=None, field=None, timesteps=None, subset=None, av
 		# Special case in 2D: we have to prepare for pcolormesh instead of imshow
 		elif len(self._centers) == 2:
 			p1 = self._centers[0] # locations of grid points along first dimension
-			d = (p1[1,:] - p1[0,:])/2. # half separation between the points
-			p1[0,:] += d
-			p1 = self._np.vstack((p1, p1[-1,:]+d)) # add last edges at the end of box
-			offset = p1[0,:]
-			p1 = self._np.apply_along_axis(lambda x: x-offset, 1, p1) # move points
 			p2 = self._centers[1] # locations of grid points along second dimension
-			d = (p2[1,:] - p2[0,:])/2. # half separation between the points
-			p2[0,:] += d
-			p2 = self._np.vstack((p2, p2[-1,:]+d)) # add last edges at the end of box
-			offset = p2[0,:]
-			p2 = self._np.apply_along_axis(lambda x: x-offset, 1, p2) # move points
+			d1 = (p1[1,:] - p1[0,:]) # separation between the points
+			d2 = (p2[1,:] - p2[0,:])
+			p1 = self._np.vstack((p1, p1[-1,:]+d1)) # add last edges at the end of box
+			p2 = self._np.vstack((p2, p2[-1,:]+d2))
+			p1 -= 0.5*(d1+d2) # Move all edges by half separation
+			p2 -= 0.5*(d1+d2)
 			# Trick in a 3D simulation (the probe has to be projected)
 			if self._ndim_particles==3:
 				# unit vectors in the two dimensions + perpendicular
 				u1 = self.p_plot[0] / self._np.linalg.norm(self.p_plot[0])
 				u2 = self.p_plot[1] / self._np.linalg.norm(self.p_plot[1])
+				# Prepare offset (zero = box origin, projected on the probe plane)
+				u1u2 = self._np.dot(u1, u2)
+				Ox = self._np.dot(p1[0,:], u1)
+				Oy = (self._np.dot(p1[0,:], u2)-Ox*u1u2)/(1-u1u2**2)
 				# Distances along first direction
-				p1[:,0] = self._np.dot(p1, u1)
-				p1[:,1:] = 0.
+				p1[:,0] = self._np.dot(p1-p1[0,:], u1) + Ox
+				p1[:,1] = Oy
+				p1[:,2] = 0.
 				# Distances along second direction
-				p2x = self._np.dot(p2, u1)
-				p2[:,1] = self._np.dot(p2, u2)
-				p2[:,0] = p2x
+				p2x = self._np.dot(p2-p2[0,:], u1)
+				p2[:,1] = self._np.dot(p2-p2[0,:], u2) + Oy
+				p2[:,0] = p2x + Ox
 				p2[:,2:] = 0.
 			# Now p1 and p2 contain edges grid points along the 2 dimensions
 			# We have to convert into X and Y 2D arrays (similar to meshgrid)
@@ -193,11 +194,11 @@ def _init(self, requestedProbe=None, field=None, timesteps=None, subset=None, av
 			for i in range(p2.shape[0]):
 				X[:,i] = p1[:,0] + (p2[i,0]-p2[0,0])
 				Y[:,i] = p1[:,1] + (p2[i,1]-p2[0,1])
-			#FOLLOWING LINES CREATE PB IN THE SELECTION OF min/max OF AXES
-			#X = self._np.maximum( X, 0.)
-			#X = self._np.minimum( X, self._ncels[0]*self._cell_length[0])
-			#Y = self._np.maximum( Y, 0.)
-			#Y = self._np.minimum( Y, self._ncels[1]*self._cell_length[1])
+			if self._ndim_particles==2:
+				X = self._np.maximum( X, 0.)
+				X = self._np.minimum( X, self._ncels[0]*self._cell_length[0])
+				Y = self._np.maximum( Y, 0.)
+				Y = self._np.minimum( Y, self._ncels[1]*self._cell_length[1])
 			self._edges = [X, Y]
 			self._limits = [[X.min(), X.max()],[Y.min(), Y.max()]]
 
@@ -346,12 +347,14 @@ def limits(self):
 		"""
 		assert self.dim <= 2, "Method limits() may only be used in 1D or 2D"
 		self._prepare1()
-		l = []
 		factor = [self._xfactor, self._yfactor]
-		for i in range(self.dim):
-			l.append([self._limits[i][0]*factor[i], self._limits[i][1]*factor[i]])
+		offset = [self._xoffset, self._yoffset]
+		l = [[
+				(offset[i] + self._limits[i][0])*factor[i], 
+				(offset[i] + self._limits[i][1])*factor[i]
+			] for i in range(self.dim)]
 		return l
-	
+
 	# get all available fields
 	def getFields(self):
 		return self._fields
@@ -414,7 +417,10 @@ def _prepare4(self):
 
 	# Overloading a plotting function in order to use pcolormesh instead of imshow
 	def _plotOnAxes_2D_(self, ax, A):
-		self._plot = ax.pcolormesh(self._xfactor*self._edges[0], self._yfactor*self._edges[1], A, **self.options.image)
+		self._plot = ax.pcolormesh(
+			self._xfactor*( self._xoffset + self._edges[0] ),
+			self._yfactor*( self._yoffset + self._edges[1] ),
+			A, **self.options.image)
 		return self._plot
 	def _animateOnAxes_2D_(self, ax, A):
 		self._plot.set_array( A.flatten() )