Merge branch 'develop' of into develop

doc/Sphinx/Overview/material.rst

@@ -46,6 +46,13 @@ As of November 2021, 90 papers have been published covering a broad range of top
 .. READ THIS 
    There is now a utility to add new entries to this list.
    Use the python script doc/doi2publications.py to generate entries from a DOI number, and paste them here
+
+.. [Jonnerby2023]
+
+      J. Jonnerby, A. von Boetticher, J. Holloway, L. Corner, A. Picksley, A. J. Ross, R. J. Shalloo , C. Thornton, N. Bourgeois, R. Walczak, and S. M. Hooker,
+      `Measurement of the decay of laser-driven linear plasma wakefields`,
+      `Phys. Rev. E 108, 055211  (2023) <https://link.aps.org/doi/10.1103/PhysRevE.108.055211>`_
+         
 .. [Drobniak2023]
 
       P. Drobniak, E. Baynard, C. Bruni, K. Cassou, C. Guyot, G. Kane, S. Kazamias, V. Kubytskyi, N. Lericheux, B. Lucas, M. Pittman, F. Massimo, A. Beck, A. Specka, P. Nghiem, and D. Minenna,

doc/Sphinx/Understand/azimuthal_modes_decomposition.rst

@@ -384,13 +384,13 @@ Equations :eq:`transverse_on_axis` can have a non zero solution only for :math:`
 We therefore conclude that all modes must cancel on axis except for :math:`m=1`.
 
 .. math::
-   E_{\theta}^{m>1}[2] = 0
+   \tilde{E_{\theta}}^{m>1}[2] = 0
 
-   B_r^{m>1}[2] = 0
+   \tilde{B_r}^{m>1}[2] = 0
 
-   E_r^{m>1}[2] = -E_r^{m>1}[3]
+   \tilde{E_r}^{m>1}[2] = -\tilde{E_r}^{m>1}[3]
 
-   B_{\theta}^{m>1}[2] = -B_{\theta}^{m>1}[3]
+   \tilde{B_{\theta}}^{m>1}[2] = -\tilde{B_{\theta}}^{m>1}[3]
 
 Let's now write the Gauss law for mode :math:`m=1`:
 
@@ -407,49 +407,19 @@ The continuity equation on axis and written in cylindrical coordinates becomes:
 
 Eq. :eq:`transverse_on_axis` already establishes that the first term is zero.
 It is only necessary to cancel the second term.
-
-In order to do so, let's build an uncentered finite dfference scheme of the second order.
-Simple Taylor developments give for any quantity :math:`u`:
-
-.. math::
-   u(x+\frac{dx}{2})=u(x)+\frac{dx}{2}u'(x)+\frac{dx^2}{8}u''(x)+O(dx3)
-
-   u(x+\frac{3dx}{2})=u(x)+\frac{3dx}{2}u'(x)+\frac{9dx^2}{8}u''(x)+O(dx3)
-
-By combination we obtain the scheme we are looking for:
-
-.. math::
-   u'(x) = \frac{9u(x+\frac{dx}{2})-u(x+\frac{3dx}{2})-8u(x)}{3dx}
-
-We can therefore write:
-
-.. math::
-   \frac{\partial \tilde{E_r}^{m=1}}{\partial r}(r=0)= 9\tilde{E_r}^{m=1}(r=\frac{dr}{2})-\tilde{E_r}^{m=1}(r=\frac{3dr}{2})-8\tilde{E_r}^{m=1}(r=0) = 0
-
-which gives:
-
-.. math::
-   \tilde{E_r}^{m=1}(r=0)=\frac{1}{8}\left(9\tilde{E_r}^{m=1}(r=\frac{dr}{2})-\tilde{E_r}^{m=1}(r=\frac{3dr}{2})\right)
-
-And from :eq:`transverse_on_axis` this turns into:
-
-.. math::
-   \tilde{E_{\theta}}^{m=1}(r=0)=\frac{-i}{8}\left(9\tilde{E_r}^{m=1}(r=\frac{dr}{2})-\tilde{E_r}^{m=1}(r=\frac{3dr}{2})\right)
-
-giving the corresponding boundary condition for :math:`E_{\theta}^{m=1}`:
+To ensure this derivative cancels on axis we simply pick:
 
 .. math::
-   E_{\theta}^{m=1}[2] = \frac{-i}{8}\left(9E_r^{m=1}[3]-E_r^{m=1}[4]\right)
+   \tilde{E_r}^{m=1}[2] = \tilde{E}_r^{m=1}[3]
 
-Once :math:`E_{\theta}^{m=1}` is defined on axis,  we need to pick :math:`E_r^{m=1}` so that :eq:`transverse_on_axis` is matched.
-With a linear interpolation we obtain:
+And equation :eq:`transverse_on_axis` then gives 
 
 .. math::
-   E_r^{m=1}[2] = 2iE_{\theta}^{m=1}[2]-E_r^{m=1}[3]
+   \tilde{E_{\theta}}^{m=1}[2] = -i\tilde{E_r}^{m=1}[3]
 
 All the equation derived here are also valid for the magnetic field.
 But because of a different duality, it is more convenient to use a different approach.
-The equations :eq:`MaxwellEqsAzimuthalModes` has a :math:`\frac{E_l}{r}` term in the expression of :math:`B_r` which makes it undefined on axis.
+The equations :eq:`MaxwellEqsAzimuthalModes` has a :math:`\frac{\tilde{E_l}}{r}` term in the expression of :math:`B_r` which makes it undefined on axis.
 Nevertheless, we need to evaluate this term for the mode :math:`m=1` and it can be done as follows.
 
 .. math::

happi/_Diagnostics/Diagnostic.py

@@ -129,6 +129,7 @@ def limits(self):
 		--------
 		A list of [min, max] for each axis.
 		"""
+		assert self.dim <= 2, "Method limits() may only be used in 1D or 2D"
 		self._prepare1()
 		l = []
 		factor = [self._xfactor, self._yfactor]

happi/_Diagnostics/Probe.py

@@ -154,9 +154,11 @@ def _init(self, requestedProbe=None, field=None, timesteps=None, subset=None, av
 
 		self._selection = tuple(s if type(s) is slice else slice(s,s+1) for s in self._selection)
 
+
 		# Special case in 1D: we convert the point locations to scalar distances
 		if len(self._centers) == 1:
 			self._centers[0] = self._np.sqrt(self._np.sum((self._centers[0]-self._centers[0][0])**2,axis=1))
+			self._limits = [[self._centers[0].min(), self._centers[0].max()]]
 		# 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
@@ -197,6 +199,7 @@ def _init(self, requestedProbe=None, field=None, timesteps=None, subset=None, av
 			#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()]]
 
 		# Prepare the reordering of the points for patches disorder
 		tmpShape = self._initialShape
@@ -333,6 +336,22 @@ def changeField(self, field):
 		self._loadField(field)
 		self._prepareUnits()
 
+	# Method to obtain the plot limits
+	def limits(self):
+		"""Gets the overall limits of the diagnostic along its axes
+
+		Returns:
+		--------
+		A list of [min, max] for each axis.
+		"""
+		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]])
+		return l
+
 	# get all available fields
 	def getFields(self):
 		return self._fields

src/ElectroMagnSolver/MA_SolverAM_norm.cpp

@@ -34,8 +34,6 @@ void MA_SolverAM_norm::operator()( ElectroMagn *fields )
         cField2D *Jt = ( static_cast<ElectroMagnAM *>( fields ) )->Jt_[imode];
         int j_glob    = ( static_cast<ElectroMagnAM *>( fields ) )->j_glob_;
         bool isYmin = ( static_cast<ElectroMagnAM *>( fields ) )->isYmin;
-        //double *invR = ( static_cast<ElectroMagnAM *>( fields ) )->invR;
-        //double *invRd = ( static_cast<ElectroMagnAM *>( fields ) )->invRd;
 
         // Electric field Elr^(d,p)
         for( unsigned int i=0 ; i<nl_d ; i++ ) {
@@ -69,7 +67,6 @@ void MA_SolverAM_norm::operator()( ElectroMagn *fields )
                     ( *Et )( i, j-1 )=-( *Et )( i, j+1 );
                 }
                 for( unsigned int i=0 ; i<nl_p  ; i++ ) {
-                    //( *Er )( i, j+1 )= ( *Er )( i, j+2 ) / 9.;
                     ( *Er )( i, j )= -( *Er )( i, j+1 );
                 }
                 for( unsigned int i=0 ; i<nl_d ; i++ ) {
@@ -82,20 +79,18 @@ void MA_SolverAM_norm::operator()( ElectroMagn *fields )
                     ( *El )( i, j-1 )=-( *El )( i, j+1 );
                 }
                 for( unsigned int i=0 ; i<nl_p  ; i++ ) {
-                    //( *Et )( i, j )= -1./3.*( 4.*Icpx*( *Er )( i, j+1 )+( *Et )( i, j+1 ) );
-                    ( *Et )( i, j )= -Icpx/8.*( 9.*( *Er )( i, j+1 )-( *Er )( i, j+2 ) );
+                    ( *Et )( i, j )= -Icpx*( *Er )( i, j+1 );
                     ( *Et )( i, j-1 )=( *Et )( i, j+1 );
                 }
                 for( unsigned int i=0 ; i<nl_p ; i++ ) {
-                    ( *Er )( i, j )=2.*Icpx*( *Et )( i, j )-( *Er )( i, j+1 );
+                    ( *Er )( i, j ) = ( *Er )( i, j+1 );
                 }
             } else { // mode > 1
                 for( unsigned int  i=0 ; i<nl_d; i++ ) {
                     ( *El )( i, j )= 0;
                     ( *El )( i, j-1 )=-( *El )( i, j+1 );
                 }
                 for( unsigned int  i=0 ; i<nl_p; i++ ) {
-                    ( *Er )( i, j+1 )= ( *Er )( i, j+2 ) / 9.;
                     ( *Er )( i, j )= -( *Er )( i, j+1 );
                 }
                 for( unsigned int i=0 ; i<nl_p; i++ ) {

src/ElectroMagnSolver/MF_SolverAM_Yee.cpp

@@ -44,8 +44,6 @@ void MF_SolverAM_Yee::operator()( ElectroMagn *fields )
         cField2D *Bt = ( static_cast<ElectroMagnAM *>( fields ) )->Bt_[imode];
         int  j_glob = ( static_cast<ElectroMagnAM *>( fields ) )->j_glob_;
         bool isYmin = ( static_cast<ElectroMagnAM *>( fields ) )->isYmin;
-        //double *invR = ( static_cast<ElectroMagnAM *>( fields ) )->invR;
-        //double *invRd = ( static_cast<ElectroMagnAM *>( fields ) )->invRd;
 
         // Magnetic field Bl^(p,d)
         for( unsigned int i=0 ; i<nl_p;  i++ ) {
@@ -81,7 +79,6 @@ void MF_SolverAM_Yee::operator()( ElectroMagn *fields )
                     ( *Br )( i, 1 )=-( *Br )( i, 3 );
                 }
                 for( unsigned int i=0 ; i<nl_d ; i++ ) {
-                    //( *Bt )( i, j+1 )= ( *Bt )( i, j+2 )/9.;
                     ( *Bt )( i, j )= -( *Bt )( i, j+1 );
                 }
                 for( unsigned int i=0 ; i<nl_p ; i++ ) {
@@ -100,7 +97,7 @@ void MF_SolverAM_Yee::operator()( ElectroMagn *fields )
                     ( *Br )( i, 1 )=( *Br )( i, 3 );
                 }
                 for( unsigned int i=0; i<nl_d ; i++ ) {
-                    ( *Bt )( i, j )= -2.*Icpx*( *Br )( i, j )-( *Bt )( i, j+1 );
+                    ( *Bt )( i, j )= ( *Bt )( i, j+1 );
                 }
 
             } else { // modes > 1

src/Field/Field3D.cpp

@@ -274,7 +274,7 @@ double Field3D::norm2OnDevice( unsigned int istart[3][2], unsigned int bufsize[3
 #if defined( SMILEI_ACCELERATOR_GPU_OMP )
     #pragma omp target teams distribute parallel for collapse(3) \
 	      map(tofrom: nrj)  \
-              map(from: ixstart, ixend, iystart, iyend, izstart, izend) \
+              map(to: ny, nz, ixstart, ixend, iystart, iyend, izstart, izend) \
 	      /*is_device_ptr( data_ ) */           \
 	      reduction(+:nrj) 
 #elif defined( SMILEI_OPENACC_MODE )

src/Projector/ProjectorAM2Order.cpp

@@ -460,22 +460,23 @@ void ProjectorAM2Order::axisBC(ElectroMagnAM *emAM, bool diag_flag )
 
 void ProjectorAM2Order::apply_axisBC(std::complex<double> *rhoj,std::complex<double> *Jl, std::complex<double> *Jr, std::complex<double> *Jt, unsigned int imode, bool diag_flag )
 {
-
-   double sign = -1.;
+   double sign = 1.;
    for (unsigned int i=0; i< imode; i++) sign *= -1;
 
    if (diag_flag && rhoj) {
        for( unsigned int i=2 ; i<npriml_*nprimr_+2; i+=nprimr_ ) {
            //Fold rho
            for( unsigned int j=1 ; j<3; j++ ) {
                rhoj[i+j] += sign * rhoj[i-j];
-               rhoj[i-j]  = sign * rhoj[i+j];
            }
            //Apply BC
            if (imode > 0){
                rhoj[i] = 0.;
+               rhoj[i-1]  = - rhoj[i+1];
            } else {
+	       //This is just for cosmetics on the picture, rho has no influence on the results
                rhoj[i] = (4.*rhoj[i+1] - rhoj[i+2])/3.;
+               rhoj[i-1]  = rhoj[i+1];
            }
        }
    }
@@ -484,14 +485,14 @@ void ProjectorAM2Order::apply_axisBC(std::complex<double> *rhoj,std::complex<dou
        for( unsigned int i=2 ; i<(npriml_+1)*nprimr_+2; i+=nprimr_ ) {
            //Fold Jl
            for( unsigned int j=1 ; j<3; j++ ) {
-               Jl [i+j] +=  sign * Jl[i-j];
-               Jl[i-j]   =  sign * Jl[i+j];
+               Jl [i+j] +=  sign * Jl[i-j]; //Add even modes, substract odd modes since el(theta=0 = el(theta=pi) at all r.
             }
             if (imode > 0){
                 Jl [i] = 0. ;
+                Jl[i-1]   =  -Jl[i+1];
            } else {
-                //Force dJl/dr = 0 at r=0.
-                Jl [i] =  (4.*Jl [i+1] - Jl [i+2])/3. ;
+		//Jl mode 1 on axis is left as is. It looks over estimated but it is necessary to conserve a correct divergence and a proper evaluation on the field on axis.
+                Jl [i-1] =  Jl [i+1] ;
            }
        }
    }
@@ -502,24 +503,19 @@ void ProjectorAM2Order::apply_axisBC(std::complex<double> *rhoj,std::complex<dou
            int ilocr = i*(nprimr_+1)+3;
            //Fold Jt
            for( unsigned int j=1 ; j<3; j++ ) {
-               Jt [iloc+j] += -sign * Jt[iloc-j];
-               Jt[iloc-j]   = -sign * Jt[iloc+j];
+               Jt [iloc+j] -= sign * Jt[iloc-j]; // substract pair modes, add odd modes since et(theta=0 = -et(theta=pi) at all r.
            }
            for( unsigned int j=0 ; j<3; j++ ) {
-               Jr [ilocr+2-j] += -sign * Jr [ilocr-3+j];
-               Jr[ilocr-3+j]     = -sign * Jr[ilocr+2-j];
+               Jr [ilocr+2-j] -= sign * Jr [ilocr-3+j];// substract pair modes, add odd modes since er(theta=0 = -er(theta=pi) at all r.
            }
 
            if (imode == 1){
-               Jt [iloc]= -Icpx/8.*( 9.*Jr[ilocr]- Jr[ilocr+1]);
-               //Force dJr/dr = 0 at r=0.
-               //Jr [ilocr] =  (25.*Jr[ilocr+1] - 9*Jr[ilocr+2])/16. ;
-               Jr [ilocr-1] = 2.*Icpx*Jt[iloc] - Jr [ilocr];
+               Jr [ilocr-1] = Jr [ilocr]; // dJr/dr mode 1 = 0 on axis
+               Jt [iloc]= -Icpx*Jr[ilocr]; // Jt mode 1 = -I Jr mode 1 on axis to keep div(J) = 0.
            } else{
-               Jt [iloc] = 0. ;
-               //Force dJr/dr = 0 and Jr=0 at r=0.
-               //Jr [ilocr] =  Jr [ilocr+1]/9.;
-               Jr [ilocr-1] = -Jr [ilocr];
+               Jt [iloc] = 0. ; // only mode 1 is non zero on axis
+               Jt [iloc-1] = -Jt [iloc+1]; // only mode 1 is non zero on axis
+               Jr [ilocr-1] = -Jr [ilocr]; // only mode 1 is non zero on axis
            }
        }
    }

validation/easi/__init__.py

@@ -556,6 +556,7 @@ def list_benchmarks(self):
 
         return benchmarks
 
+    #Compare the results of given "simulation_dir" to the reference results of the given benchname.
     def compare(self, benchname, simulation_dir):
         from os import chdir
         from .tools import execfile
@@ -569,6 +570,18 @@ def compare(self, benchname, simulation_dir):
             print("The validation procedure failed.")
         else:
             print("PASS")
+
+    #Generate reference from an already existing directory
+    def generate_reference(self, benchname, simulation_dir):
+        from os import chdir
+        from .tools import execfile
+        global _dataNotMatching
+        _dataNotMatching = False
+        validation_script = self.smilei_path.analyses + "validate_" + benchname
+        Validate = self.CreateReference(self.smilei_path.references, benchname)
+        chdir(simulation_dir)
+        execfile(validation_script, {"Validate":Validate})
+        Validate.write()
 
     # DEFINE A CLASS TO CREATE A REFERENCE
     class CreateReference(object):