Merge branch 'v8.2.0'

002_change_version_number.py

@@ -9,7 +9,7 @@
         with open(filename) as fid:
             content = fid.read()
         if 'giovanni.iadarola@cern.ch' in content:
-            content = content.replace('PyECLOUD Version 8.1.0', 'PyECLOUD Version 8.1.0')
+            content = content.replace('PyECLOUD Version 8.2.0', 'PyECLOUD Version 8.2.0')
             with open(filename, 'w') as fid:
                 fid.write(content)
 

003_change_preamble.py

@@ -20,7 +20,7 @@
 #
 #     This file is part of the code:
 #
-#                   PyECLOUD Version 8.1.0
+#                   PyECLOUD Version 8.2.0
 #
 #
 #     Main author:          Giovanni IADAROLA
@@ -32,6 +32,7 @@
 #
 #     Contributors:         Eleonora Belli
 #                           Philipp Dijkstal
+#                           Lorenzo Giacomel
 #                           Lotta Mether
 #                           Annalisa Romano
 #                           Giovanni Rumolo

MP_system.py

@@ -7,7 +7,7 @@
 #
 #     This file is part of the code:
 #
-#                   PyECLOUD Version 8.1.0
+#                   PyECLOUD Version 8.2.0
 #
 #
 #     Main author:          Giovanni IADAROLA
@@ -19,6 +19,7 @@
 #
 #     Contributors:         Eleonora Belli
 #                           Philipp Dijkstal
+#                           Lorenzo Giacomel
 #                           Lotta Mether
 #                           Annalisa Romano
 #                           Giovanni Rumolo
@@ -68,7 +69,7 @@ def __init__(self, N_mp_max, nel_mp_ref_0, fact_split, fact_clean,
                  Dx_hist_reg, Nx_reg, Ny_reg, Nvx_reg, Nvy_reg, Nvz_reg, regen_hist_cut, chamb,
                  N_mp_soft_regen=None, N_mp_after_soft_regen=None,
                  N_mp_async_regen=None, N_mp_after_async_regen=None,
-                 charge=-e, mass=m_e, flag_lifetime_hist = False):
+                 charge=-e, mass=m_e, flag_lifetime_hist = False, name=None):
 
         N_mp_max = int(N_mp_max)
         self.x_mp = np.zeros(N_mp_max, float)
@@ -107,6 +108,7 @@ def __init__(self, N_mp_max, nel_mp_ref_0, fact_split, fact_clean,
 
         self.charge = charge
         self.mass = mass
+        self.name = name
 
         xg_hist_reg = np.arange(0, chamb.x_aper + 2. * Dx_hist_reg, Dx_hist_reg, float)
         xgr_hist_reg = xg_hist_reg[1:]
@@ -132,7 +134,7 @@ def __init__(self, N_mp_max, nel_mp_ref_0, fact_split, fact_clean,
 
     def clean_small_MPs(self):
 
-        print "Start clean. N_mp=%d Nel=%e"%(self.N_mp, np.sum(self.nel_mp[0:self.N_mp]))
+        print "Cloud %s: Start clean. N_mp=%d Nel=%e"%(self.name, self.N_mp, np.sum(self.nel_mp[0:self.N_mp]))
 
         flag_clean = (self.nel_mp < self.nel_mp_cl_th)
         flag_keep = ~(flag_clean)
@@ -152,7 +154,11 @@ def clean_small_MPs(self):
         if self.flag_lifetime_hist:
             self.t_last_impact[0:self.N_mp] = np.array(self.t_last_impact[flag_keep].copy())
 
-        print "Done clean. N_mp=%d Nel=%e"%(self.N_mp, np.sum(self.nel_mp[0:self.N_mp]))
+        print "Cloud %s: Done clean. N_mp=%d Nel=%e"%(self.name, self.N_mp, np.sum(self.nel_mp[0:self.N_mp]))
+
+        if self.N_mp == 0:
+            self.set_nel_mp_ref(self.nel_mp_ref_0)
+            print('Cloud %s: nel_mp_ref set to nel_mp_ref_0'%self.name)
 
     def set_nel_mp_ref(self, val):
         self.nel_mp_ref = val
@@ -170,7 +176,7 @@ def perform_soft_regeneration(self, target_N_mp):
                     new_nel_mp_ref = self.nel_mp_ref_0
 
                 #if new_nel_mp_ref>self.nel_mp_ref_0:removed from version 3.16
-                print 'Start SOFT regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.N_mp, chrg, erg)
+                print 'Cloud %s: Start SOFT regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.name, self.N_mp, chrg, erg)
 
                 self.set_nel_mp_ref(new_nel_mp_ref)
 
@@ -198,13 +204,13 @@ def perform_soft_regeneration(self, target_N_mp):
 
                 correct_fact = chrg_before / chrg_after
 
-                print 'Applied correction factor = %e'%correct_fact
+                print 'Cloud %s: Applied correction factor = %e'%(self.name, correct_fact)
 
                 self.nel_mp[0:self.N_mp] = self.nel_mp[0:self.N_mp] * correct_fact
 
                 chrg = np.sum(self.nel_mp)
                 erg = np.sum(0.5 / np.abs(self.charge / self.mass) * self.nel_mp[0:self.N_mp] * (self.vx_mp[0:self.N_mp] * self.vx_mp[0:self.N_mp] + self.vy_mp[0:self.N_mp] * self.vy_mp[0:self.N_mp] + self.vz_mp[0:self.N_mp] * self.vz_mp[0:self.N_mp]))
-                print 'Done SOFT regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.N_mp, chrg, erg)
+                print 'Cloud %s: Done SOFT regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.name, self.N_mp, chrg, erg)
 
     def check_for_soft_regeneration(self):
 
@@ -216,15 +222,17 @@ def check_for_async_regeneration(self):
 
         if self.flag_async_regen:
             if self.N_mp > self.N_mp_async_regen:
-                print('Asynchronous regeneration.')
-                self.perform_soft_regeneration(target_N_mp=self.N_mp_after_async_regen)
+                print('Cloud %s: Asynchronous clean and regeneration.' %self.name)
+                self.clean_small_MPs()
+                if self.N_mp > self.N_mp_async_regen:
+                    self.perform_soft_regeneration(target_N_mp=self.N_mp_after_async_regen)
 
     def check_for_regeneration(self):
 
         if (self.N_mp > self.N_mp_regen or (self.N_mp < self.N_mp_regen_low and self.nel_mp_ref > self.nel_mp_ref_0)):
             chrg = np.sum(self.nel_mp)
             erg = np.sum(0.5 / np.abs(self.charge / self.mass) * self.nel_mp[0:self.N_mp] * (self.vx_mp[0:self.N_mp] * self.vx_mp[0:self.N_mp] + self.vy_mp[0:self.N_mp] * self.vy_mp[0:self.N_mp] + self.vz_mp[0:self.N_mp] * self.vz_mp[0:self.N_mp]))
-            print 'Start regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.N_mp, chrg, erg)
+            print 'Cloud %s: Start regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.name, self.N_mp, chrg, erg)
 
             new_nel_mp_ref = chrg / self.N_mp_after_regen
             if new_nel_mp_ref < self.nel_mp_ref_0:
@@ -245,7 +253,7 @@ def check_for_regeneration(self):
 
             x_max = (len(hist_vect) - i_cut + 1) * self.Dx_hist_reg + self.bias_x_hist_reg
 
-            print 'x_max = %e'%x_max
+            print 'Cloud %s: x_max = %e'%(self.name, x_max)
 
             flag_clean = (abs(self.x_mp) > x_max)
             flag_keep = ~(flag_clean)
@@ -295,7 +303,7 @@ def check_for_regeneration(self):
             bias_vz = np.ceil(float(self.Nvz_reg) / 2)
             #Attention when trnslating to python
 
-            print 'particles_assigned_to grid'
+            print 'Cloud %s: particles_assigned_to grid'%(self.name)
 
             ##
             #% MATLAB-like indices
@@ -428,7 +436,7 @@ def check_for_regeneration(self):
 
             chrg = np.sum(self.nel_mp)
             erg = np.sum(0.5 / np.abs(self.charge / self.mass) * self.nel_mp[0:self.N_mp] * (self.vx_mp[0:self.N_mp] * self.vx_mp[0:self.N_mp] + self.vy_mp[0:self.N_mp] * self.vy_mp[0:self.N_mp] + self.vz_mp[0:self.N_mp] * self.vz_mp[0:self.N_mp]))
-            print 'Done regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.N_mp, chrg, erg)
+            print 'Cloud %s: Done regeneration. N_mp=%d Nel_tot=%1.2e En_tot=%1.2e'%(self.name, self.N_mp, chrg, erg)
 
     def add_uniform_MP_distrib(self, DNel, E_init, x_max, x_min, y_max, y_min):
 

PyEC4PyHT.py

@@ -7,7 +7,7 @@
 #
 #     This file is part of the code:
 #
-#                   PyECLOUD Version 8.1.0
+#                   PyECLOUD Version 8.2.0
 #
 #
 #     Main author:          Giovanni IADAROLA
@@ -19,6 +19,7 @@
 #
 #     Contributors:         Eleonora Belli
 #                           Philipp Dijkstal
+#                           Lorenzo Giacomel
 #                           Lotta Mether
 #                           Annalisa Romano
 #                           Giovanni Rumolo
@@ -104,7 +105,7 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean
                  beam_monitor=None, verbose=False, save_pyecl_outp_as=None,
                  **kwargs):
 
-        print 'PyECLOUD Version 8.1.0'
+        print 'PyECLOUD Version 8.2.0'
 
         # These git commands return the hash and the branch of the specified git directory.
         path_to_git = os.path.dirname(os.path.abspath(__file__)) + '/.git'

beam_and_timing.py

@@ -7,7 +7,7 @@
 #
 #     This file is part of the code:
 #
-#                   PyECLOUD Version 8.1.0
+#                   PyECLOUD Version 8.2.0
 #
 #
 #     Main author:          Giovanni IADAROLA
@@ -19,6 +19,7 @@
 #
 #     Contributors:         Eleonora Belli
 #                           Philipp Dijkstal
+#                           Lorenzo Giacomel
 #                           Lotta Mether
 #                           Annalisa Romano
 #                           Giovanni Rumolo

boris_c_function.c

@@ -4,15 +4,18 @@ void boris_c(int N_sub_steps, double Dtt,
 		double* B_field, double* B_skew,
 		double* xn1, double* yn1,  double* zn1,
 		double* vxn1, double* vyn1, double* vzn1,
-		double* Ex_n, double* Ey_n, int N_mp, int N_multipoles,
+		double* Ex_n, double* Ey_n, 
+		double* Bx_n_custom, double* By_n_custom, double* Bz_n_custom,
+		int custom_B,
+		int N_mp, int N_multipoles,
 		double charge, double mass)
 {
 	int p, isub, order;
 	double Ex_np, Ey_np;
-	double Bx_n, By_n;
+	double Bx_n, By_n, Bz_n;
 	double rexy, imxy, rexy_0;
-	double tBx, tBy, tBsq;
-	double sBx, sBy;
+	double tBx, tBy, tBz, tBsq;
+	double sBx, sBy, sBz;
 	double vx_prime, vy_prime, vz_prime;
 	double vx_min, vy_min, vz_min;
 	double vx_plus, vy_plus, vz_plus;
@@ -39,6 +42,7 @@ void boris_c(int N_sub_steps, double Dtt,
 			imxy = 0.;
 			By_n = B_field[0];
 			Bx_n = B_skew[0];
+			Bz_n = 0.;
 
 			for(order = 1; order < N_multipoles; order++)
 			{
@@ -54,27 +58,35 @@ void boris_c(int N_sub_steps, double Dtt,
 				By_n += (B_field[order]*rexy - B_skew[order]*imxy);
 				Bx_n += (B_field[order]*imxy + B_skew[order]*rexy);
 			}
+
+			if (custom_B){
+				Bx_n += Bx_n_custom[p];
+				By_n += By_n_custom[p]; 
+				Bz_n += Bz_n_custom[p]; 
+			}
 
 			tBx = 0.5*qm*Dtt*Bx_n;
 			tBy = 0.5*qm*Dtt*By_n;
-			tBsq = tBx*tBx + tBy*tBy;
+			tBz = 0.5*qm*Dtt*Bz_n;
+			tBsq = tBx*tBx + tBy*tBy + tBz*tBz;
 
 			sBx = 2.*tBx/(1.+tBsq);
 			sBy = 2.*tBy/(1.+tBsq);
+			sBz = 2.*tBz/(1.+tBsq);
 
 			vx_min = vxn1p + 0.5*qm*Ex_np*Dtt;
 			vy_min = vyn1p + 0.5*qm*Ey_np*Dtt;
 			vz_min = vzn1p;
 
 			//v_prime = v_min + cross(v_min, tB)
-			vx_prime = -vz_min*tBy + vx_min;
-			vy_prime = vz_min*tBx + vy_min;
-			vz_prime = vx_min*tBy-vy_min*tBx + vz_min;
+			vx_prime = vy_min*tBz - vz_min*tBy + vx_min;
+			vy_prime = vz_min*tBx - vx_min*tBz + vy_min;
+			vz_prime = vx_min*tBy - vy_min*tBx + vz_min;
 
 			//v_plus = v_min + cross(v_prime, sB)
-			vx_plus = -vz_prime*sBy + vx_min;
-			vy_plus = vz_prime*sBx + vy_min;
-			vz_plus = vx_prime*sBy-vy_prime*sBx + vz_min;
+			vx_plus = vy_prime*sBz - vz_prime*sBy + vx_min;
+			vy_plus = vz_prime*sBx - vx_prime*sBz+ vy_min;
+			vz_plus = vx_prime*sBy - vy_prime*sBx + vz_min;
 
 			vxn1p = vx_plus + 0.5*qm*Ex_np*Dtt;
 			vyn1p = vy_plus + 0.5*qm*Ey_np*Dtt;

boris_c_function.h

@@ -5,7 +5,10 @@ void boris_c(int N_sub_steps, double Dtt,
 		double* B_field, double* B_skew,
 		double* xn1, double* yn1,  double* zn1, 
 		double* vxn1, double* vyn1, double* vzn1,
-		double* Ex_n, double* Ey_n, int N_mp, int N_multipoles,
+		double* Ex_n, double* Ey_n,
+		double* Bx_n_custom, double* By_n_custom, double* Bz_n_custom,
+		int custom_B,
+		int N_mp, int N_multipoles,
 		double charge, double mass);
 
 #endif

boris_cython.pxd

@@ -3,6 +3,8 @@ cdef extern from "boris_c_function.h":
     void boris_c(int N_sub_steps, double Dtt, double* B_field, double* B_skew,
                           double* xn1, double* yn1,  double* zn1,
                           double* vxn1, double* vyn1, double* vzn1,
-                          double* Ex_n, double* Ey_n, int N_mp, int N_multipoles,
+                          double* Ex_n, double* Ey_n,
+                          double* Bx_n_custom, double* By_n_custom, double* Bz_n_custom,
+		                      int custom_B,
+                          int N_mp, int N_multipoles,
                           double charge, double mass)
-

boris_cython.pyx

@@ -6,7 +6,10 @@ cimport numpy as np
 cpdef boris_step_multipole(N_sub_steps, Dtt, np.ndarray B_field, np.ndarray B_skew,
                            np.ndarray xn1, np.ndarray yn1, np.ndarray zn1,
                            np.ndarray vxn1, np.ndarray vyn1, np.ndarray vzn1,
-                           np.ndarray Ex_n, np.ndarray Ey_n, charge, mass):
+                           np.ndarray Ex_n, np.ndarray Ey_n,
+                           np.ndarray Bx_n, np.ndarray By_n, np.ndarray Bz_n,
+                           custom_B,
+                           charge, mass):
 
 
     cdef double* B_field_data = <double*>B_field.data
@@ -19,10 +22,14 @@ cpdef boris_step_multipole(N_sub_steps, Dtt, np.ndarray B_field, np.ndarray B_sk
     cdef double* vzn1_data =  <double*>vzn1.data
     cdef double* Ex_n_data =  <double*>Ex_n.data
     cdef double* Ey_n_data =  <double*>Ey_n.data
+    cdef double* Bx_n_data =  <double*>Bx_n.data
+    cdef double* By_n_data =  <double*>By_n.data
+    cdef double* Bz_n_data =  <double*>Bz_n.data
 
     boris_c(N_sub_steps, Dtt, B_field_data, B_skew_data,
           xn1_data, yn1_data,  zn1_data,
           vxn1_data, vyn1_data, vzn1_data,
-          Ex_n_data, Ey_n_data, len(xn1), len(B_field),
+          Ex_n_data, Ey_n_data,
+          Bx_n_data, By_n_data, Bz_n_data, custom_B,
+          len(xn1), len(B_field),
           charge, mass)
-