diff --git a/002_change_version_number.py b/002_change_version_number.py index b29cf671..8139c96e 100644 --- a/002_change_version_number.py +++ b/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.2.0', 'PyECLOUD Version 8.2.0') + content = content.replace('PyECLOUD Version 8.3.0', 'PyECLOUD Version 8.3.0') with open(filename, 'w') as fid: fid.write(content) diff --git a/003_change_preamble.py b/003_change_preamble.py index 71a4eaab..564b5d71 100644 --- a/003_change_preamble.py +++ b/003_change_preamble.py @@ -20,7 +20,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/MP_system.py b/MP_system.py index 6e7d55b3..624911d9 100644 --- a/MP_system.py +++ b/MP_system.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/PyEC4PyHT.py b/PyEC4PyHT.py index 080d1616..724bd313 100644 --- a/PyEC4PyHT.py +++ b/PyEC4PyHT.py @@ -1,4 +1,4 @@ -#-Begin-preamble------------------------------------------------------- +# -Begin-preamble------------------------------------------------------- # # CERN # @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA @@ -48,17 +48,16 @@ # The material cannot be sold. CERN should be given credit in # all references. # -#-End-preamble--------------------------------------------------------- +# -End-preamble--------------------------------------------------------- import os import subprocess import time import numpy as np -from scipy.constants import c, e, m_e +from scipy.constants import c from . import myloadmat_to_obj as mlm -from . import init from . import buildup_simulation as bsim @@ -70,63 +69,100 @@ class DummyBeamTim(object): def __init__(self, PyPIC_state): self.PyPIC_state = PyPIC_state - self.b_spac = 0. + self.b_spac = 0.0 self.pass_numb = 0 self.N_pass_tot = 1 def get_beam_eletric_field(self, MP_e): - if (MP_e.N_mp > 0): + if MP_e.N_mp > 0: if self.PyPIC_state is None: - Ex_n_beam = 0. * MP_e.x_mp[0:MP_e.N_mp] - Ey_n_beam = 0. * MP_e.y_mp[0:MP_e.N_mp] + Ex_n_beam = 0.0 * MP_e.x_mp[0 : MP_e.N_mp] + Ey_n_beam = 0.0 * MP_e.y_mp[0 : MP_e.N_mp] else: # compute beam electric field - Ex_n_beam, Ey_n_beam = self.PyPIC_state.gather(MP_e.x_mp[0:MP_e.N_mp], MP_e.y_mp[0:MP_e.N_mp]) + Ex_n_beam, Ey_n_beam = self.PyPIC_state.gather( + MP_e.x_mp[0 : MP_e.N_mp], MP_e.y_mp[0 : MP_e.N_mp] + ) else: - Ex_n_beam = 0. - Ey_n_beam = 0. + Ex_n_beam = 0.0 + Ey_n_beam = 0.0 return Ex_n_beam, Ey_n_beam -extra_allowed_kwargs = {'x_beam_offset', 'y_beam_offset', 'probes_position', 'enable_kick_x', 'enable_kick_y'} +extra_allowed_kwargs = { + "x_beam_offset", + "y_beam_offset", + "probes_position", + "enable_kick_x", + "enable_kick_y", +} class Ecloud(object): - def generate_twin_ecloud_with_shared_space_charge(self): - if hasattr(self, 'efieldmap'): - raise ValueError('Ecloud has been replaced with field map. I cannot generate a twin ecloud!') - return Ecloud(self.L_ecloud, self.slicer, self.Dt_ref, self.pyecl_input_folder, self.flag_clean_slices, - self.slice_by_slice_mode, self.spacech_ele, self.kick_mode_for_beam_field, - self.beam_monitor, self.verbose, self.save_pyecl_outp_as, **self.kwargs) - - def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean_slices=False, - slice_by_slice_mode=False, space_charge_obj=None, kick_mode_for_beam_field=False, - beam_monitor=None, verbose=False, save_pyecl_outp_as=None, - **kwargs): - - print('PyECLOUD Version 8.2.0') + if hasattr(self, "efieldmap"): + raise ValueError( + "Ecloud has been replaced with field map. I cannot generate a twin ecloud!" + ) + return Ecloud( + self.L_ecloud, + self.slicer, + self.Dt_ref, + self.pyecl_input_folder, + self.flag_clean_slices, + self.slice_by_slice_mode, + self.spacech_ele, + self.kick_mode_for_beam_field, + self.beam_monitor, + self.verbose, + self.save_pyecl_outp_as, + self.force_interp_at_substeps_interacting_slices, + **self.kwargs + ) + + def __init__( + self, + L_ecloud, + slicer, + Dt_ref, + pyecl_input_folder="./", + flag_clean_slices=False, + slice_by_slice_mode=False, + space_charge_obj=None, + kick_mode_for_beam_field=False, + beam_monitor=None, + verbose=False, + save_pyecl_outp_as=None, + force_interp_at_substeps_interacting_slices=False, + **kwargs + ): + + print("PyECLOUD Version 8.3.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' - cmd_hash = 'git --git-dir %s rev-parse HEAD' % path_to_git - cmd_branch = 'git --git-dir %s rev-parse --abbrev-ref HEAD' % path_to_git + path_to_git = os.path.dirname(os.path.abspath(__file__)) + "/.git" + cmd_hash = "git --git-dir %s rev-parse HEAD" % path_to_git + cmd_branch = "git --git-dir %s rev-parse --abbrev-ref HEAD" % path_to_git try: - git_hash = 'git hash: %s' % (subprocess.check_output(cmd_hash.split()).split()[0]) + git_hash = "git hash: %s" % ( + subprocess.check_output(cmd_hash.split()).split()[0] + ) except Exception as e: - git_hash = 'Retrieving git hash failed' + git_hash = "Retrieving git hash failed" print(e) print(git_hash) try: - git_branch = 'git branch: %s' % (subprocess.check_output(cmd_branch.split()).split()[0]) + git_branch = "git branch: %s" % ( + subprocess.check_output(cmd_branch.split()).split()[0] + ) except Exception as e: - git_branch = 'Retrieving git branch failed' + git_branch = "Retrieving git branch failed" print(e) print(git_branch) - print('PyHEADTAIL module') - print('Initializing ecloud from folder: ' + pyecl_input_folder) + print("PyHEADTAIL module") + print("Initializing ecloud from folder: " + pyecl_input_folder) self.slicer = slicer self.Dt_ref = Dt_ref self.L_ecloud = L_ecloud @@ -134,38 +170,48 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean self.pyecl_input_folder = pyecl_input_folder self.kwargs = kwargs + self.force_interp_at_substeps_interacting_slices = ( + force_interp_at_substeps_interacting_slices + ) + + print(f"Reinterp fields at substeps for interacting slices: {self.force_interp_at_substeps_interacting_slices}") + self.cloudsim = bsim.BuildupSimulation( - pyecl_input_folder=pyecl_input_folder, skip_beam=True, + pyecl_input_folder=pyecl_input_folder, + skip_beam=True, spacech_ele=space_charge_obj, ignore_kwargs=extra_allowed_kwargs, skip_pyeclsaver=(save_pyecl_outp_as is None), filen_main_outp=save_pyecl_outp_as, - **self.kwargs) + **self.kwargs + ) - if self.cloudsim.config_dict['track_method'] == 'Boris': + if self.cloudsim.config_dict["track_method"] == "Boris": pass - elif self.cloudsim.config_dict['track_method'] == 'BorisMultipole': + elif self.cloudsim.config_dict["track_method"] == "BorisMultipole": pass else: - raise ValueError("""track_method should be 'Boris' or 'BorisMultipole' - others are not implemented in the PyEC4PyHT module""") + raise ValueError( + """track_method should be 'Boris' or 'BorisMultipole' - others are not implemented in the PyEC4PyHT module""" + ) - self.x_beam_offset = 0. - self.y_beam_offset = 0. - if 'x_beam_offset' in kwargs: - self.x_beam_offset = kwargs['x_beam_offset'] - if 'y_beam_offset' in kwargs: - self.y_beam_offset = kwargs['y_beam_offset'] + self.x_beam_offset = 0.0 + self.y_beam_offset = 0.0 + if "x_beam_offset" in kwargs: + self.x_beam_offset = kwargs["x_beam_offset"] + if "y_beam_offset" in kwargs: + self.y_beam_offset = kwargs["y_beam_offset"] self.enable_kick_x = True self.enable_kick_y = True - if 'enable_kick_x' in kwargs: - self.enable_kick_x = kwargs['enable_kick_x'] + if "enable_kick_x" in kwargs: + self.enable_kick_x = kwargs["enable_kick_x"] if not self.enable_kick_x: - print('Horizontal kick on the beam is disabled!') - if 'enable_kick_y' in kwargs: - self.enable_kick_y = kwargs['enable_kick_y'] + print("Horizontal kick on the beam is disabled!") + if "enable_kick_y" in kwargs: + self.enable_kick_y = kwargs["enable_kick_y"] if not self.enable_kick_y: - print('Vertical kick on the beam is disabled!') + print("Vertical kick on the beam is disabled!") # initialize proton density probes self.save_ele_field_probes = False @@ -173,25 +219,25 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean self.y_probes = -1 self.Ex_ele_last_track_at_probes = -1 self.Ey_ele_last_track_at_probes = -1 - if 'probes_position' in list(kwargs.keys()): + if "probes_position" in list(kwargs.keys()): self.save_ele_field_probes = True - self.probes_position = kwargs['probes_position'] + self.probes_position = kwargs["probes_position"] self.N_probes = len(self.probes_position) self.x_probes = [] self.y_probes = [] for ii_probe in range(self.N_probes): - self.x_probes.append(self.probes_position[ii_probe]['x']) - self.y_probes.append(self.probes_position[ii_probe]['y']) + self.x_probes.append(self.probes_position[ii_probe]["x"]) + self.y_probes.append(self.probes_position[ii_probe]["y"]) self.x_probes = np.array(self.x_probes) self.y_probes = np.array(self.y_probes) self.N_tracks = 0 - #self.cloudsim.spacech_ele.flag_decimate = False + # self.cloudsim.spacech_ele.flag_decimate = False if self.cloudsim.flag_multiple_clouds: - raise ValueError('Multiple clouds not yet tested in PyEC4PyHT!') + raise ValueError("Multiple clouds not yet tested in PyEC4PyHT!") self.save_ele_distributions_last_track = False self.save_ele_potential_and_field = False @@ -200,7 +246,7 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean self.save_ele_MP_position = False self.save_ele_MP_velocity = False self.save_ele_MP_size = False - + self.save_beam_distributions_last_track = False self.save_beam_potential_and_field = False self.save_beam_potential = False @@ -208,7 +254,9 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean self.track_only_first_time = False - self.initial_MP_e_clouds = [cl.MP_e.extract_dict() for cl in self.cloudsim.cloud_list] + self.initial_MP_e_clouds = [ + cl.MP_e.extract_dict() for cl in self.cloudsim.cloud_list + ] self.flag_clean_slices = flag_clean_slices @@ -219,7 +267,7 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean self.track = self._track_in_single_slice_mode self.finalize_and_reinitialize = self._finalize_and_reinitialize - self.spacech_ele = self.cloudsim.spacech_ele # For backwards compatibility + self.spacech_ele = self.cloudsim.spacech_ele # For backwards compatibility self.kick_mode_for_beam_field = kick_mode_for_beam_field self.beam_monitor = beam_monitor @@ -228,7 +276,7 @@ def __init__(self, L_ecloud, slicer, Dt_ref, pyecl_input_folder='./', flag_clean self.save_pyecl_outp_as = save_pyecl_outp_as self.i_reinit = 0 - self.t_sim = 0. + self.t_sim = 0.0 self.i_curr_bunch = -1 # @profile @@ -236,7 +284,7 @@ def track(self, beam): if self.track_only_first_time: if self.N_tracks > 0: - print('Warning: Track skipped because track_only_first_time is True.') + print("Warning: Track skipped because track_only_first_time is True.") return if self.verbose: @@ -244,8 +292,8 @@ def track(self, beam): self._reinitialize() - if hasattr(beam.particlenumber_per_mp, '__iter__'): - raise ValueError('ecloud module assumes same size for all beam MPs') + if hasattr(beam.particlenumber_per_mp, "__iter__"): + raise ValueError("ecloud module assumes same size for all beam MPs") if self.flag_clean_slices: beam.clean_slices() @@ -254,13 +302,13 @@ def track(self, beam): for i in range(slices.n_slices - 1, -1, -1): if self.verbose: - print(('Slice %d/%d'%(i, slices.n_slices))) + print(("Slice %d/%d" % (i, slices.n_slices))) # select particles in the slice ix = slices.particle_indices_of_slice(i) # slice size and time step - dz = (slices.z_bins[i + 1] - slices.z_bins[i]) + dz = slices.z_bins[i + 1] - slices.z_bins[i] self._track_single_slice(beam, ix, dz, force_pyecl_newpass=(i == 0)) @@ -272,24 +320,33 @@ def track(self, beam): if self.verbose: stop_time = time.mktime(time.localtime()) - print('Done track %d in %.1f s'%(self.N_tracks, stop_time - start_time)) + print("Done track %d in %.1f s" % (self.N_tracks, stop_time - start_time)) self.N_tracks += 1 def replace_with_recorded_field_map(self, delete_ecloud_data=True): if self.track_only_first_time: - print('Warning: replace_with_recorded_field_map resets track_only_first_time = False') + print( + "Warning: replace_with_recorded_field_map resets track_only_first_time = False" + ) self.track_only_first_time = False - if not hasattr(self, 'efieldmap'): + if not hasattr(self, "efieldmap"): from .Transverse_Efield_map_for_frozen_cloud import Transverse_Efield_map - self.efieldmap = Transverse_Efield_map(xg=self.spacech_ele.xg, yg=self.spacech_ele.yg, - Ex=self.Ex_ele_last_track, Ey=self.Ey_ele_last_track, L_interaction=self.L_ecloud, - slicer=self.slicer, - flag_clean_slices=True, - x_beam_offset=self.x_beam_offset, y_beam_offset=self.y_beam_offset, - slice_by_slice_mode=self.slice_by_slice_mode) + + self.efieldmap = Transverse_Efield_map( + xg=self.spacech_ele.xg, + yg=self.spacech_ele.yg, + Ex=self.Ex_ele_last_track, + Ey=self.Ey_ele_last_track, + L_interaction=self.L_ecloud, + slicer=self.slicer, + flag_clean_slices=True, + x_beam_offset=self.x_beam_offset, + y_beam_offset=self.y_beam_offset, + slice_by_slice_mode=self.slice_by_slice_mode, + ) self._ecloud_track = self.track @@ -301,14 +358,16 @@ def replace_with_recorded_field_map(self, delete_ecloud_data=True): self.cloudsim = None else: - print('Warning: efieldmap already exists. I do nothing.') + print("Warning: efieldmap already exists. I do nothing.") - def track_once_and_replace_with_recorded_field_map(self, bunch, delete_ecloud_data=True): + def track_once_and_replace_with_recorded_field_map( + self, bunch, delete_ecloud_data=True + ): self.save_ele_field = True self.track_only_first_time = True if self.slice_by_slice_mode: - if not hasattr(bunch, '__iter__'): - raise ValueError('A list of slices should be provided!') + if not hasattr(bunch, "__iter__"): + raise ValueError("A list of slices should be provided!") self._reinitialize() for slc in bunch: self.track(slc) @@ -324,9 +383,9 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): spacech_ele = self.cloudsim.spacech_ele # Check if the slice interacts with the beam - if hasattr(slic, 'slice_info'): - if 'interact_with_EC' in list(slic.slice_info.keys()): - interact_with_EC = slic.slice_info['interact_with_EC'] + if hasattr(slic, "slice_info"): + if "interact_with_EC" in list(slic.slice_info.keys()): + interact_with_EC = slic.slice_info["interact_with_EC"] else: interact_with_EC = True else: @@ -336,12 +395,16 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): dt_slice = dz / (slic.beta * c) # Check if sub-slicing is needed - if self.cloudsim.config_dict['Dt'] is not None: - if dt_slice > self.cloudsim.config_dict['Dt']: + if self.cloudsim.config_dict["Dt"] is not None: + if dt_slice > self.cloudsim.config_dict["Dt"]: if interact_with_EC: - raise ValueError('Slices that interact with the cloud cannot be longer than the buildup timestep!') + raise ValueError( + "Slices that interact with the cloud cannot be longer than the buildup timestep!" + ) - N_cloud_steps = np.int_(np.ceil(dt_slice / self.cloudsim.config_dict['Dt'])) + N_cloud_steps = np.int_( + np.ceil(dt_slice / self.cloudsim.config_dict["Dt"]) + ) dt_cloud_step = dt_slice / N_cloud_steps dt_array = np.array(N_cloud_steps * [dt_cloud_step]) else: @@ -350,16 +413,19 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): dt_array = np.array([dt_slice]) # Acquire bunch passage information - if hasattr(slic, 'slice_info'): - if 'info_parent_bunch' in list(slic.slice_info.keys()): + if hasattr(slic, "slice_info"): + if "info_parent_bunch" in list(slic.slice_info.keys()): # check if first slice of first bunch - if slic.slice_info['info_parent_bunch']['i_bunch'] == 0 and slic.slice_info['i_slice'] == 0: + if ( + slic.slice_info["info_parent_bunch"]["i_bunch"] == 0 + and slic.slice_info["i_slice"] == 0 + ): self.finalize_and_reinitialize() # check if new passage - if slic.slice_info['info_parent_bunch']['i_bunch'] > self.i_curr_bunch: - self.i_curr_bunch = slic.slice_info['info_parent_bunch']['i_bunch'] + if slic.slice_info["info_parent_bunch"]["i_bunch"] > self.i_curr_bunch: + self.i_curr_bunch = slic.slice_info["info_parent_bunch"]["i_bunch"] new_pass = True else: new_pass = force_pyecl_newpass @@ -381,44 +447,49 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): N_sub_steps = 1 Dt_substep = dt / N_sub_steps - #print Dt_substep, N_sub_steps, dt + # print Dt_substep, N_sub_steps, dt if len(ix) == 0 or not interact_with_EC: # no particles in the beam - #build dummy beamtim object + # build dummy beamtim object dummybeamtim = DummyBeamTim(None) - dummybeamtim.lam_t_curr = 0. - dummybeamtim.sigmax = 0. - dummybeamtim.sigmay = 0. - dummybeamtim.x_beam_pos = 0. - dummybeamtim.y_beam_pos = 0. + dummybeamtim.lam_t_curr = 0.0 + dummybeamtim.sigmax = 0.0 + dummybeamtim.sigmay = 0.0 + dummybeamtim.x_beam_pos = 0.0 + dummybeamtim.y_beam_pos = 0.0 else: # beam field self.beam_PyPIC_state.scatter( x_mp=slic.x[ix] + self.x_beam_offset, y_mp=slic.y[ix] + self.y_beam_offset, - nel_mp=slic.x[ix] * 0. + slic.particlenumber_per_mp / dz, - charge=slic.charge) + nel_mp=slic.x[ix] * 0.0 + slic.particlenumber_per_mp / dz, + charge=slic.charge, + ) self.cloudsim.spacech_ele.PyPICobj.solve_states([self.beam_PyPIC_state]) - #build dummy beamtim object + # build dummy beamtim object dummybeamtim = DummyBeamTim(self.beam_PyPIC_state) - dummybeamtim.lam_t_curr = np.mean(slic.particlenumber_per_mp / dz) * len(ix) + dummybeamtim.lam_t_curr = np.mean( + slic.particlenumber_per_mp / dz + ) * len(ix) dummybeamtim.sigmax = np.std(slic.x[ix]) dummybeamtim.sigmay = np.std(slic.y[ix]) dummybeamtim.x_beam_pos = np.mean(slic.x[ix]) + self.x_beam_offset dummybeamtim.y_beam_pos = np.mean(slic.y[ix]) + self.y_beam_offset - dummybeamtim.tt_curr = self.t_sim # In order to have the PIC activated + dummybeamtim.tt_curr = self.t_sim # In order to have the PIC activated dummybeamtim.Dt_curr = dt dummybeamtim.pass_numb = self.i_curr_bunch dummybeamtim.flag_new_bunch_pass = new_pass # Force space charge recomputation - force_recompute_space_charge = interact_with_EC or (i_clou_step == 0) # we always force at the first step, as we don't know the state of the PIC + force_recompute_space_charge = interact_with_EC or ( + i_clou_step == 0 + ) # we always force at the first step, as we don't know the state of the PIC # Disable cleanings and regenerations skip_MP_cleaning = interact_with_EC @@ -427,11 +498,17 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): # print(dummybeamtim.tt_curr, dummybeamtim.flag_new_bunch_pass, force_recompute_space_charge) # Perform cloud simulation step - self.cloudsim.sim_time_step(beamtim_obj=dummybeamtim, - Dt_substep_custom=Dt_substep, N_sub_steps_custom=N_sub_steps, - kick_mode_for_beam_field=self.kick_mode_for_beam_field, - force_recompute_space_charge=force_recompute_space_charge, - skip_MP_cleaning=skip_MP_cleaning, skip_MP_regen=skip_MP_regen) + self.cloudsim.sim_time_step( + beamtim_obj=dummybeamtim, + Dt_substep_custom=Dt_substep, + N_sub_steps_custom=N_sub_steps, + kick_mode_for_beam_field=self.kick_mode_for_beam_field, + force_recompute_space_charge=force_recompute_space_charge, + skip_MP_cleaning=skip_MP_cleaning, + skip_MP_regen=skip_MP_regen, + force_reinterp_fields_at_substeps=(interact_with_EC + and self.force_interp_at_substeps_interacting_slices) + ) if interact_with_EC: # Build MP_system-like object with beam coordinates @@ -444,7 +521,11 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): Ex_sc_p, Ey_sc_p = spacech_ele.get_sc_eletric_field(MP_p) ## kick beam particles - fact_kick = slic.charge / (slic.mass * slic.beta * slic.beta * slic.gamma * c * c) * self.L_ecloud + fact_kick = ( + slic.charge + / (slic.mass * slic.beta * slic.beta * slic.gamma * c * c) + * self.L_ecloud + ) if self.enable_kick_x: slic.xp[ix] += fact_kick * Ex_sc_p if self.enable_kick_y: @@ -484,14 +565,18 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): if self.save_ele_MP_size: self.nel_MP_last_track.append(MPe_for_save.nel_mp.copy()) - if self.save_ele_MP_position or self.save_ele_MP_velocity or self.save_ele_MP_size: + if ( + self.save_ele_MP_position + or self.save_ele_MP_velocity + or self.save_ele_MP_size + ): self.N_MP_last_track.append(MPe_for_save.N_mp) if self.save_ele_field_probes: MP_probes = Empty() MP_probes.x_mp = self.x_probes MP_probes.y_mp = self.y_probes - MP_probes.nel_mp = self.x_probes * 0. + 1. # fictitious charge of 1 C + MP_probes.nel_mp = self.x_probes * 0.0 + 1.0 # fictitious charge of 1 C MP_probes.N_mp = len(self.x_probes) Ex_sc_probe, Ey_sc_probe = spacech_ele.get_sc_eletric_field(MP_probes) @@ -499,13 +584,15 @@ def _track_single_slice(self, slic, ix, dz, force_pyecl_newpass=False): self.Ey_ele_last_track_at_probes.append(Ey_sc_probe.copy()) self.t_sim += dt - new_pass = False # it can be true only for the first sub-slice of a slice + new_pass = False # it can be true only for the first sub-slice of a slice def _reinitialize(self): cc = mlm.obj_from_dict(self.cloudsim.config_dict) - for thiscloud, initdict in zip(self.cloudsim.cloud_list, self.initial_MP_e_clouds): + for thiscloud, initdict in zip( + self.cloudsim.cloud_list, self.initial_MP_e_clouds + ): thiscloud.MP_e.init_from_dict(initdict) thiscloud.MP_e.set_nel_mp_ref(thiscloud.MP_e.nel_mp_ref_0) @@ -513,19 +600,45 @@ def _reinitialize(self): if thiscloud.pyeclsaver is not None: thiscloud.pyeclsaver.extract_sey = False - thiscloud.pyeclsaver.start_observing(cc.Dt, thiscloud.MP_e, None, thiscloud.impact_man, - thisconf.r_center, thisconf.Dt_En_hist, thisconf.logfile_path, thisconf.progress_path, - flag_detailed_MP_info=thisconf.flag_detailed_MP_info, flag_movie=thisconf.flag_movie, - flag_sc_movie=thisconf.flag_sc_movie, save_mp_state_time_file=thisconf.save_mp_state_time_file, - flag_presence_sec_beams=self.cloudsim.flag_presence_sec_beams, sec_beams_list=self.cloudsim.sec_beams_list, dec_fac_secbeam_prof=thisconf.dec_fac_secbeam_prof, - el_density_probes=thisconf.el_density_probes, save_simulation_state_time_file=thisconf.save_simulation_state_time_file, - x_min_hist_det=thisconf.x_min_hist_det, x_max_hist_det=thisconf.x_max_hist_det, - y_min_hist_det=thisconf.y_min_hist_det, y_max_hist_det=thisconf.y_max_hist_det, - Dx_hist_det=thisconf.Dx_hist_det, dec_fact_out=cc.dec_fact_out, stopfile=cc.stopfile, filen_main_outp=thisconf.filen_main_outp, - flag_cos_angle_hist=thisconf.flag_cos_angle_hist, cos_angle_width=thisconf.cos_angle_width, - flag_multiple_clouds=(len(self.cloudsim.cloud_list) > 1), cloud_name=thisconf.cloud_name, flag_last_cloud=(thiscloud is self.cloudsim.cloud_list[-1])) - - thiscloud.pyeclsaver.filen_main_outp = thiscloud.pyeclsaver.filen_main_outp.split('.mat')[0].split('__iter')[0] + '__iter%d.mat'%self.i_reinit + thiscloud.pyeclsaver.start_observing( + cc.Dt, + thiscloud.MP_e, + None, + thiscloud.impact_man, + thisconf.r_center, + thisconf.Dt_En_hist, + thisconf.logfile_path, + thisconf.progress_path, + flag_detailed_MP_info=thisconf.flag_detailed_MP_info, + flag_movie=thisconf.flag_movie, + flag_sc_movie=thisconf.flag_sc_movie, + save_mp_state_time_file=thisconf.save_mp_state_time_file, + flag_presence_sec_beams=self.cloudsim.flag_presence_sec_beams, + sec_beams_list=self.cloudsim.sec_beams_list, + dec_fac_secbeam_prof=thisconf.dec_fac_secbeam_prof, + el_density_probes=thisconf.el_density_probes, + save_simulation_state_time_file=thisconf.save_simulation_state_time_file, + x_min_hist_det=thisconf.x_min_hist_det, + x_max_hist_det=thisconf.x_max_hist_det, + y_min_hist_det=thisconf.y_min_hist_det, + y_max_hist_det=thisconf.y_max_hist_det, + Dx_hist_det=thisconf.Dx_hist_det, + dec_fact_out=cc.dec_fact_out, + stopfile=cc.stopfile, + filen_main_outp=thisconf.filen_main_outp, + flag_cos_angle_hist=thisconf.flag_cos_angle_hist, + cos_angle_width=thisconf.cos_angle_width, + flag_multiple_clouds=(len(self.cloudsim.cloud_list) > 1), + cloud_name=thisconf.cloud_name, + flag_last_cloud=(thiscloud is self.cloudsim.cloud_list[-1]), + ) + + thiscloud.pyeclsaver.filen_main_outp = ( + thiscloud.pyeclsaver.filen_main_outp.split(".mat")[0].split( + "__iter" + )[0] + + "__iter%d.mat" % self.i_reinit + ) if self.save_ele_distributions_last_track: self.rho_ele_last_track = [] @@ -566,14 +679,18 @@ def _reinitialize(self): if self.save_ele_MP_size: self.nel_MP_last_track = [] - if self.save_ele_MP_position or self.save_ele_MP_velocity or self.save_ele_MP_size: + if ( + self.save_ele_MP_position + or self.save_ele_MP_velocity + or self.save_ele_MP_size + ): self.N_MP_last_track = [] if self.save_ele_field_probes: self.Ex_ele_last_track_at_probes = [] self.Ey_ele_last_track_at_probes = [] - #~ self.t_sim = 0. + # ~ self.t_sim = 0. self.i_curr_bunch = -1 self.i_reinit += 1 @@ -588,7 +705,7 @@ def _finalize(self): if self.save_ele_field: self.Ex_ele_last_track = np.array(self.Ex_ele_last_track[::-1]) self.Ey_ele_last_track = np.array(self.Ey_ele_last_track[::-1]) - + if self.save_beam_distributions_last_track: self.rho_beam_last_track = np.array(self.rho_beam_last_track[::-1]) @@ -610,33 +727,40 @@ def _finalize(self): if self.save_ele_MP_size: self.nel_MP_last_track = np.array(self.nel_MP_last_track[::-1]) - if self.save_ele_MP_position or self.save_ele_MP_velocity or self.save_ele_MP_size: - self.N_MP_last_track = np.array(self.N_MP_last_track[::-1]) + if ( + self.save_ele_MP_position + or self.save_ele_MP_velocity + or self.save_ele_MP_size + ): + self.N_MP_last_track = np.array(self.N_MP_last_track[::-1]) if self.save_ele_field_probes: - self.Ex_ele_last_track_at_probes = np.array(self.Ex_ele_last_track_at_probes[::-1]) - self.Ey_ele_last_track_at_probes = np.array(self.Ey_ele_last_track_at_probes[::-1]) + self.Ex_ele_last_track_at_probes = np.array( + self.Ex_ele_last_track_at_probes[::-1] + ) + self.Ey_ele_last_track_at_probes = np.array( + self.Ey_ele_last_track_at_probes[::-1] + ) def _finalize_and_reinitialize(self): - print('Exec. finalize and reinitialize') + print("Exec. finalize and reinitialize") self._finalize() self._reinitialize() def _track_in_single_slice_mode(self, beam): - if hasattr(beam.particlenumber_per_mp, '__iter__'): - raise ValueError('ecloud module assumes same size for all beam MPs') + if hasattr(beam.particlenumber_per_mp, "__iter__"): + raise ValueError("ecloud module assumes same size for all beam MPs") if self.flag_clean_slices: - raise ValueError( - 'track cannot clean the slices in slice-by-slice mode! ') + raise ValueError("track cannot clean the slices in slice-by-slice mode! ") - if beam.slice_info != 'unsliced': - dz = beam.slice_info['z_bin_right'] - beam.slice_info['z_bin_left'] - self._track_single_slice(beam, ix=np.arange(beam.macroparticlenumber), dz=dz) + if beam.slice_info != "unsliced": + dz = beam.slice_info["z_bin_right"] - beam.slice_info["z_bin_left"] + self._track_single_slice( + beam, ix=np.arange(beam.macroparticlenumber), dz=dz + ) def remove_savers(self): for thiscloud in self.cloudsim.cloud_list: thiscloud.pyeclsaver = None - - diff --git a/beam_and_timing.py b/beam_and_timing.py index ba182870..faa9e5e4 100644 --- a/beam_and_timing.py +++ b/beam_and_timing.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/buildup_simulation.py b/buildup_simulation.py index b8e132dd..3ae0cad3 100644 --- a/buildup_simulation.py +++ b/buildup_simulation.py @@ -1,6 +1,6 @@ #!/afs/cern.ch/project/uslarp/opt/lxplus64/Python-2.7.2/bin/python -#-Begin-preamble------------------------------------------------------- +# -Begin-preamble------------------------------------------------------- # # CERN # @@ -9,7 +9,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA @@ -50,7 +50,7 @@ # The material cannot be sold. CERN should be given credit in # all references. # -#-End-preamble--------------------------------------------------------- +# -End-preamble--------------------------------------------------------- from . import init as init @@ -60,21 +60,39 @@ class BuildupSimulation(object): - - def __init__(self, pyecl_input_folder='./', skip_beam=False, skip_spacech_ele=False, - skip_pyeclsaver=False, ignore_kwargs=[], spacech_ele=None, **kwargs): - - print('PyECLOUD Version 8.2.0') - beamtim, spacech_ele, t_sc_ON, flag_presence_sec_beams, sec_beams_list, \ - config_dict, flag_multiple_clouds, cloud_list, checkpoint_folder, \ - cross_ion = init.read_input_files_and_init_components(\ - pyecl_input_folder=pyecl_input_folder, - skip_beam=skip_beam, - skip_pyeclsaver=skip_pyeclsaver, - skip_spacech_ele=skip_spacech_ele, - spacech_ele=spacech_ele, - ignore_kwargs=ignore_kwargs, - **kwargs) + def __init__( + self, + pyecl_input_folder="./", + skip_beam=False, + skip_spacech_ele=False, + skip_pyeclsaver=False, + ignore_kwargs=[], + spacech_ele=None, + **kwargs + ): + + print("PyECLOUD Version 8.3.0") + ( + beamtim, + spacech_ele, + t_sc_ON, + flag_presence_sec_beams, + sec_beams_list, + config_dict, + flag_multiple_clouds, + cloud_list, + checkpoint_folder, + cross_ion, + flag_reinterp_fields_at_substeps, + ) = init.read_input_files_and_init_components( + pyecl_input_folder=pyecl_input_folder, + skip_beam=skip_beam, + skip_pyeclsaver=skip_pyeclsaver, + skip_spacech_ele=skip_spacech_ele, + spacech_ele=spacech_ele, + ignore_kwargs=ignore_kwargs, + **kwargs + ) self.config_dict = config_dict self.beamtim = beamtim @@ -89,19 +107,26 @@ def __init__(self, pyecl_input_folder='./', skip_beam=False, skip_spacech_ele=Fa self.flag_em_tracking = spacech_ele.flag_em_tracking self.cross_ion = cross_ion + self.flag_reinterp_fields_at_substeps = flag_reinterp_fields_at_substeps + print(f"Reinterp fields at substeps: {self.flag_reinterp_fields_at_substeps}") + # Checking if there are saved checkpoints if self.checkpoint_folder is not None: if os.path.isdir(self.checkpoint_folder): if len(os.listdir(self.checkpoint_folder)) == 1: - print('Loading from checkpoint...') + print("Loading from checkpoint...") # Selecting latest checkpoint checkpoint = os.listdir(self.checkpoint_folder)[0] - self.load_checkpoint(filename_simulation_checkpoint=checkpoint, load_from_folder=self.checkpoint_folder) + self.load_checkpoint( + filename_simulation_checkpoint=checkpoint, + load_from_folder=self.checkpoint_folder, + ) elif len(os.listdir(self.checkpoint_folder)) == 0: - print('No checkpoint found, starting new simulation...') + print("No checkpoint found, starting new simulation...") else: - raise ValueError('More than one checkpoint found in %s'%self.checkpoint_folder) - + raise ValueError( + "More than one checkpoint found in %s" % self.checkpoint_folder + ) def run(self, t_end_sim=None): @@ -110,14 +135,14 @@ def run(self, t_end_sim=None): flag_presence_sec_beams = self.flag_presence_sec_beams sec_beams_list = self.sec_beams_list - print('Start timestep iter') + print("Start timestep iter") ## simulation while not beamtim.end_simulation(): if t_end_sim is not None and beamtim.tt_curr is not None: if beamtim.tt_curr >= t_end_sim: - print('Reached user defined t_end_sim --> Ending simulation') + print("Reached user defined t_end_sim --> Ending simulation") break beamtim.next_time_step() @@ -129,152 +154,309 @@ def run(self, t_end_sim=None): self.sim_time_step() if beamtim.flag_new_bunch_pass: - print('**** Done pass_numb = %d/%d\n'%(beamtim.pass_numb, beamtim.N_pass_tot)) - - - def sim_time_step(self, beamtim_obj=None, Dt_substep_custom=None, N_sub_steps_custom=None, kick_mode_for_beam_field=False, - force_recompute_space_charge=False, skip_MP_cleaning=False, skip_MP_regen=False): + print( + "**** Done pass_numb = %d/%d\n" + % (beamtim.pass_numb, beamtim.N_pass_tot) + ) + + def sim_time_step( + self, + beamtim_obj=None, + Dt_substep_custom=None, + N_sub_steps_custom=None, + kick_mode_for_beam_field=False, + force_recompute_space_charge=False, + force_reinterp_fields_at_substeps=False, + skip_MP_cleaning=False, + skip_MP_regen=False, + ): if beamtim_obj is not None: beamtim = beamtim_obj else: beamtim = self.beamtim - spacech_ele = self.spacech_ele - t_sc_ON = self.t_sc_ON - flag_presence_sec_beams = self.flag_presence_sec_beams - sec_beams_list = self.sec_beams_list - flag_multiple_clouds = self.flag_multiple_clouds - cloud_list = self.cloud_list - cross_ion = self.cross_ion - - flag_recompute_space_charge = spacech_ele.check_for_recomputation(t_curr=beamtim.tt_curr) + flag_recompute_space_charge = self.spacech_ele.check_for_recomputation( + t_curr=beamtim.tt_curr + ) # Loop over clouds: gather fields, move, generate new MPs - for i_cloud, cloud in enumerate(cloud_list): - MP_e = cloud.MP_e - dynamics = cloud.dynamics - impact_man = cloud.impact_man - gas_ion_flag = cloud.gas_ion_flag - resgasion = cloud.resgasion - t_ion = cloud.t_ion - photoem_flag = cloud.photoem_flag - phemiss = cloud.phemiss - - ## Compute beam electric field (main and secondary beams) - Ex_n_beam, Ey_n_beam = beamtim.get_beam_eletric_field(MP_e) - - if flag_presence_sec_beams: - for sec_beam in sec_beams_list: - Ex_n_secbeam, Ey_n_secbeam = sec_beam.get_beam_eletric_field(MP_e) - Ex_n_beam += Ex_n_secbeam - Ey_n_beam += Ey_n_secbeam - - ## Either compute electromagnetic or electrostatic fields - if self.flag_em_tracking: - Ex_sc_n, Ey_sc_n, Bx_sc_n, By_sc_n, Bz_sc_n = spacech_ele.get_sc_em_field(MP_e) - else: - ## Compute electron space charge electric field - Ex_sc_n, Ey_sc_n = spacech_ele.get_sc_eletric_field(MP_e) - Bx_sc_n = np.asarray([0.]) - By_sc_n = np.asarray([0.]) - Bz_sc_n = np.asarray([0.]) - - - if kick_mode_for_beam_field: - if Dt_substep_custom is None or N_sub_steps_custom is None: - raise ValueError("""Kick mode can be used only with custom time steps!""") - MP_e.vx_mp[:MP_e.N_mp] += Ex_n_beam * (Dt_substep_custom * N_sub_steps_custom) * MP_e.charge / MP_e.mass - MP_e.vy_mp[:MP_e.N_mp] += Ey_n_beam * (Dt_substep_custom * N_sub_steps_custom) * MP_e.charge / MP_e.mass - # Electric field for dynamics step - Ex_n = Ex_sc_n - Ey_n = Ey_sc_n - else: - # Electric field for dynamics step - Ex_n = Ex_sc_n + Ex_n_beam - Ey_n = Ey_sc_n + Ey_n_beam + for i_cloud, cloud in enumerate(self.cloud_list): ## Save position before motion step - old_pos = MP_e.get_positions() + old_pos = cloud.MP_e.get_positions() ## Motion - if Dt_substep_custom is None and N_sub_steps_custom is None and beamtim.flag_unif_Dt: - # Standard simulation mode - MP_e = dynamics.step(MP_e, Ex_n, Ey_n, Ez_n=0, Bx_n=Bx_sc_n, By_n=By_sc_n, Bz_n=Bz_sc_n) - elif Dt_substep_custom is None and N_sub_steps_custom is None and not(beamtim.flag_unif_Dt): - # Dt from non-uniform beam profile - if self.config_dict['track_method'] not in ['Boris', 'BorisMultipole']: - raise ValueError("""track_method should be 'Boris' or 'BorisMultipole' to use custom substeps!""") - Dt_substep_target = cloud.dynamics.Dt / cloud.dynamics.N_sub_steps - N_substeps_curr = np.round(beamtim.Dt_curr / Dt_substep_target) - Dt_substep_curr = beamtim.Dt_curr / N_substeps_curr - MP_e = dynamics.stepcustomDt(MP_e, Ex_n, Ey_n, Ez_n=0, Bx_n=Bx_sc_n, By_n=By_sc_n, Bz_n=Bz_sc_n, Dt_substep=Dt_substep_curr, N_sub_steps=N_substeps_curr) - else: - # Custom steps and substeps provided as arguments of sim_time_step - if self.config_dict['track_method'] not in ['Boris', 'BorisMultipole']: - raise ValueError("""track_method should be 'Boris' or 'BorisMultipole' to use custom substeps!""") - MP_e = dynamics.stepcustomDt(MP_e, Ex_n, Ey_n, Ez_n=0, Bx_n=Bx_sc_n, By_n=By_sc_n, Bz_n=Bz_sc_n, Dt_substep=Dt_substep_custom, N_sub_steps=N_sub_steps_custom) + ## (external B field, beam and cloud fields are taken + ## into account) + self._cloud_motion( + cloud, + beamtim, + Dt_substep_custom, + N_sub_steps_custom, + kick_mode_for_beam_field, + force_reinterp_fields_at_substeps, + ) ## Impacts: backtracking and secondary emission - MP_e = impact_man.backtrack_and_second_emiss(old_pos, MP_e, beamtim.tt_curr) + cloud.MP_e = cloud.impact_man.backtrack_and_second_emiss( + old_pos, cloud.MP_e, beamtim.tt_curr + ) ## Evolve SEY module (e.g. charge decay for insulators - impact_man.sey_mod.SEY_model_evol(Dt=beamtim.Dt_curr) - - ## Gas ionization (main and secondary beams) - if(beamtim.tt_curr < t_ion and gas_ion_flag == 1): - MP_e = resgasion.generate(MP_e, beamtim.lam_t_curr, beamtim.Dt_curr, beamtim.sigmax, beamtim.sigmay, - x_beam_pos=beamtim.x_beam_pos, y_beam_pos=beamtim.y_beam_pos) - if flag_presence_sec_beams: - for sec_beam in sec_beams_list: - MP_e = resgasion.generate(MP_e, sec_beam.lam_t_curr, sec_beam.Dt_curr, sec_beam.sigmax, sec_beam.sigmay, - x_beam_pos=sec_beam.x_beam_pos, y_beam_pos=sec_beam.y_beam_pos) - - ## Photoemission (main and secondary beams) - if (photoem_flag != 0): - lam_curr_phem = beamtim.lam_t_curr - if flag_presence_sec_beams: - for sec_beam in sec_beams_list: - lam_curr_phem += sec_beam.lam_t_curr - phemiss.generate(MP_e, lam_curr_phem, beamtim.Dt_curr) + cloud.impact_man.sey_mod.SEY_model_evol(Dt=beamtim.Dt_curr) + ## Beam-gas ionization and photoemission + self._primary_generation(cloud, beamtim) ## Cross_ionization - if cross_ion is not None: - cross_ion.generate(Dt=beamtim.Dt_curr, cloud_list=cloud_list) + if self.cross_ion is not None: + self.cross_ion.generate(Dt=beamtim.Dt_curr, cloud_list=self.cloud_list) + + ## Compute space charge field (PIC: scatter and solve) + self._recompute_cloud_spacecharge( + beamtim, flag_recompute_space_charge, force_recompute_space_charge + ) + ## Saving output + self._save_output_data(beamtim) - ## Compute space charge field - for i_cloud, cloud in enumerate(cloud_list): - if ((beamtim.tt_curr > t_sc_ON) and flag_recompute_space_charge) or force_recompute_space_charge: - flag_reset = cloud is cloud_list[0] # The first cloud resets the distribution - flag_solve = cloud is cloud_list[-1] # The last cloud computes the fields + ## Cleaning and regeneration + self._MP_cleaning_and_regenerations(beamtim, skip_MP_cleaning, skip_MP_regen) + + def _get_field_from_beams_at_particles(self, MP_e, beamtim): + Ex_n_beam, Ey_n_beam = beamtim.get_beam_eletric_field(MP_e) + + if self.flag_presence_sec_beams: + for sec_beam in self.sec_beams_list: + Ex_n_secbeam, Ey_n_secbeam = sec_beam.get_beam_eletric_field(MP_e) + Ex_n_beam += Ex_n_secbeam + Ey_n_beam += Ey_n_secbeam + return Ex_n_beam, Ey_n_beam + + def _get_field_from_clouds_at_particles(self, MP_e): + ## Either compute electromagnetic or electrostatic fields + if self.flag_em_tracking: + ( + Ex_sc_n, + Ey_sc_n, + Bx_sc_n, + By_sc_n, + Bz_sc_n, + ) = self.spacech_ele.get_sc_em_field(MP_e) + else: + ## Compute electron space charge electric field + Ex_sc_n, Ey_sc_n = self.spacech_ele.get_sc_eletric_field(MP_e) + Bx_sc_n = np.asarray([0.0]) + By_sc_n = np.asarray([0.0]) + Bz_sc_n = np.asarray([0.0]) + + return Ex_sc_n, Ey_sc_n, Bx_sc_n, By_sc_n, Bz_sc_n + + def _apply_instantaneous_kick(self, MP_e, Ex_n_kick, Ey_n_kick, Dt_kick): + + MP_e.vx_mp[: MP_e.N_mp] += Ex_n_kick * Dt_kick * MP_e.charge / MP_e.mass + MP_e.vy_mp[: MP_e.N_mp] += Ey_n_kick * Dt_kick * MP_e.charge / MP_e.mass + + def _recompute_cloud_spacecharge( + self, beamtim, flag_recompute_space_charge, force_recompute_space_charge + ): + + for i_cloud, cloud in enumerate(self.cloud_list): + if ( + (beamtim.tt_curr > self.t_sc_ON) and flag_recompute_space_charge + ) or force_recompute_space_charge: + flag_reset = ( + cloud is self.cloud_list[0] + ) # The first cloud resets the distribution + flag_solve = ( + cloud is self.cloud_list[-1] + ) # The last cloud computes the fields ## Either compute electromagnetic field or electrostatic if self.flag_em_tracking: - spacech_ele.recompute_spchg_emfield(cloud.MP_e, flag_solve=flag_solve, flag_reset=flag_reset) + self.spacech_ele.recompute_spchg_emfield( + cloud.MP_e, flag_solve=flag_solve, flag_reset=flag_reset + ) else: - spacech_ele.recompute_spchg_efield(cloud.MP_e, flag_solve=flag_solve, flag_reset=flag_reset) + self.spacech_ele.recompute_spchg_efield( + cloud.MP_e, flag_solve=flag_solve, flag_reset=flag_reset + ) # Copy rho to cloud - cloud.rho = spacech_ele.rho - sum([cl.rho for cl in cloud_list[:i_cloud]]) - - ## Saving output + cloud.rho = self.spacech_ele.rho - sum( + [cl.rho for cl in self.cloud_list[:i_cloud]] + ) + + def _cloud_motion( + self, + cloud, + beamtim, + Dt_substep_custom, + N_sub_steps_custom, + kick_mode_for_beam_field, + force_reinterp_fields_at_substeps, + ): + + flag_substeps = False + N_substeps_curr = 1 + Dt_substep_curr = None + + ## Determine mode + if N_sub_steps_custom is not None: + if self.config_dict["track_method"] not in ["Boris", "BorisMultipole"]: + raise ValueError( + """track_method should be 'Boris' or 'BorisMultipole' to use custom substeps!""" + ) + # Substepping specified as an argument of sim_time_step + flag_substeps = True + N_substeps_curr = N_sub_steps_custom + Dt_substep_curr = Dt_substep_custom + + elif not beamtim.flag_unif_Dt: + if self.config_dict["track_method"] not in ["Boris", "BorisMultipole"]: + raise ValueError( + """track_method should be 'Boris' or 'BorisMultipole' to use non-uniform timestep!""" + ) + # Non-uniform beam-profile + flag_substeps = True + Dt_substep_target = cloud.dynamics.Dt / cloud.dynamics.N_sub_steps + N_substeps_curr = int(np.round(beamtim.Dt_curr / Dt_substep_target)) + Dt_substep_curr = beamtim.Dt_curr / N_substeps_curr + + elif hasattr(cloud.dynamics, "N_sub_steps"): + # Non-unif time step is specified in the dynamics object + if cloud.dynamics.N_sub_steps is not None: + flag_substeps = True + N_substeps_curr = cloud.dynamics.N_sub_steps + Dt_substep_curr = cloud.dynamics.Dt / N_substeps_curr + + # Beam kick applided here if kick-mode (mainly used in fast-ion mode) + if kick_mode_for_beam_field: + if N_sub_steps_custom is None: + raise ValueError( + """Kick mode can be used only with custom time steps!""" + ) + + Ex_n_beam, Ey_n_beam = self._get_field_from_beams_at_particles( + cloud.MP_e, beamtim + ) + + self._apply_instantaneous_kick( + cloud.MP_e, + Ex_n_beam, + Ey_n_beam, + Dt_kick=Dt_substep_custom * N_sub_steps_custom, + ) + + # Decide number of substeps internal and external + if self.flag_reinterp_fields_at_substeps or force_reinterp_fields_at_substeps: + if not flag_substeps: + raise ValueError("No substeps set!") + N_substeps_external = N_substeps_curr + N_substeps_internal = 1 + else: + N_substeps_external = 1 + N_substeps_internal = N_substeps_curr + + # print(f'external {N_substeps_external}') + # print(f'internal {N_substeps_internal}') + for i_substep in range(N_substeps_external): + ## Interpolate fields from clouds at particles + (Ex_n, Ey_n, Bx_n, By_n, Bz_n,) = self._get_field_from_clouds_at_particles( + cloud.MP_e + ) + + ## Interpolate field from beam + if not kick_mode_for_beam_field: + Ex_n_beam, Ey_n_beam = self._get_field_from_beams_at_particles( + cloud.MP_e, beamtim + ) + Ex_n += Ex_n_beam + Ey_n += Ey_n_beam + + if not flag_substeps: + # Standard simulation mode + cloud.MP_e = cloud.dynamics.step( + cloud.MP_e, Ex_n, Ey_n, Ez_n=0, Bx_n=Bx_n, By_n=By_n, Bz_n=Bz_n, + ) + else: + # Substep mode + cloud.MP_e = cloud.dynamics.stepcustomDt( + cloud.MP_e, + Ex_n, + Ey_n, + Ez_n=0, + Bx_n=Bx_n, + By_n=By_n, + Bz_n=Bz_n, + Dt_substep=Dt_substep_curr, + N_sub_steps=N_substeps_internal, + ) + + def _primary_generation(self, cloud, beamtim): + + ## Gas ionization (main and secondary beams) + if beamtim.tt_curr < cloud.t_ion and cloud.gas_ion_flag == 1: + cloud.MP_e = cloud.resgasion.generate( + cloud.MP_e, + beamtim.lam_t_curr, + beamtim.Dt_curr, + beamtim.sigmax, + beamtim.sigmay, + x_beam_pos=beamtim.x_beam_pos, + y_beam_pos=beamtim.y_beam_pos, + ) + if self.flag_presence_sec_beams: + for sec_beam in self.sec_beams_list: + cloud.MP_e = cloud.resgasion.generate( + cloud.MP_e, + sec_beam.lam_t_curr, + sec_beam.Dt_curr, + sec_beam.sigmax, + sec_beam.sigmay, + x_beam_pos=sec_beam.x_beam_pos, + y_beam_pos=sec_beam.y_beam_pos, + ) + + ## Photoemission (main and secondary beams) + if cloud.photoem_flag != 0: + lam_curr_phem = beamtim.lam_t_curr + if self.flag_presence_sec_beams: + for sec_beam in self.sec_beams_list: + lam_curr_phem += sec_beam.lam_t_curr + cloud.phemiss.generate(cloud.MP_e, lam_curr_phem, beamtim.Dt_curr) + + def _save_output_data(self, beamtim): # We want to save and clean MP only after iteration on all clouds is completed # (e.g. to have consistent space charge state) - for cloud in cloud_list: + for cloud in self.cloud_list: if cloud.pyeclsaver is not None: # if Dt_substep_custom is not None or N_sub_steps_custom is not None: # raise ValueError('Saving with custom steps not implemented!') - cloud.impact_man = cloud.pyeclsaver.witness(cloud.MP_e, - beamtim, spacech_ele, cloud.impact_man, cloud.dynamics, - cloud.gas_ion_flag, cloud.resgasion, cloud.t_ion, t_sc_ON, - cloud.photoem_flag, cloud.phemiss, flag_presence_sec_beams, - sec_beams_list, cloud_list, buildup_sim=self, - cross_ion=self.cross_ion, rho_cloud=cloud.rho) + cloud.impact_man = cloud.pyeclsaver.witness( + cloud.MP_e, + beamtim, + self.spacech_ele, + cloud.impact_man, + cloud.dynamics, + cloud.gas_ion_flag, + cloud.resgasion, + cloud.t_ion, + self.t_sc_ON, + cloud.photoem_flag, + cloud.phemiss, + self.flag_presence_sec_beams, + self.sec_beams_list, + self.cloud_list, + buildup_sim=self, + cross_ion=self.cross_ion, + rho_cloud=cloud.rho, + ) + + def _MP_cleaning_and_regenerations(self, beamtim, skip_MP_cleaning, skip_MP_regen): - ## Cleaning and regeneration - for cloud in cloud_list: + for cloud in self.cloud_list: ## Every bunch passage if beamtim.flag_new_bunch_pass: @@ -291,23 +473,30 @@ def sim_time_step(self, beamtim_obj=None, Dt_substep_custom=None, N_sub_steps_cu cloud.MP_e.check_for_async_regeneration() + def load_state( + self, + filename_simulation_state, + force_disable_save_simulation_state=True, + filen_main_outp="Pyecltest_restarted", + load_from_folder="./", + ): # , reset_pyeclsaver = True): - def load_state(self, filename_simulation_state, force_disable_save_simulation_state=True, filen_main_outp='Pyecltest_restarted', load_from_folder='./'): # , reset_pyeclsaver = True): - - with open(load_from_folder + filename_simulation_state, 'rb') as fid: + with open(load_from_folder + filename_simulation_state, "rb") as fid: dict_state = pickle.load(fid) - self.beamtim = dict_state['beamtim'] - self.spacech_ele = dict_state['spacech_ele'] - self.t_sc_ON = dict_state['t_sc_ON'] - self.flag_presence_sec_beams = dict_state['flag_presence_sec_beams'] - self.sec_beams_list = dict_state['sec_beams_list'] + self.beamtim = dict_state["beamtim"] + self.spacech_ele = dict_state["spacech_ele"] + self.t_sc_ON = dict_state["t_sc_ON"] + self.flag_presence_sec_beams = dict_state["flag_presence_sec_beams"] + self.sec_beams_list = dict_state["sec_beams_list"] - self.flag_multiple_clouds = dict_state['flag_multiple_clouds'] + self.flag_multiple_clouds = dict_state["flag_multiple_clouds"] for i_cloud, new_cloud in enumerate(self.cloud_list): new_pyeclsaver = new_cloud.pyeclsaver - self.cloud_list[i_cloud] = dict_state['cloud_list'][i_cloud] # Replace new_cloud with saved cloud + self.cloud_list[i_cloud] = dict_state["cloud_list"][ + i_cloud + ] # Replace new_cloud with saved cloud cloud = self.cloud_list[i_cloud] # if reset_pyeclsaver or cloud.pyeclsaver is None: @@ -317,30 +506,42 @@ def load_state(self, filename_simulation_state, force_disable_save_simulation_st cloud.pyeclsaver.flag_save_simulation_state = False if filen_main_outp is not None: - filen_outp_ext = filen_main_outp.split('Pyecltest')[-1] - filen_outp_root = cloud.pyeclsaver.filen_main_outp.split('.mat')[0] - cloud.pyeclsaver.filen_main_outp = filen_outp_root + filen_outp_ext + '.mat' + filen_outp_ext = filen_main_outp.split("Pyecltest")[-1] + filen_outp_root = cloud.pyeclsaver.filen_main_outp.split(".mat")[0] + cloud.pyeclsaver.filen_main_outp = ( + filen_outp_root + filen_outp_ext + ".mat" + ) - print('Restoring PyPIC LU object...') + print("Restoring PyPIC LU object...") self.spacech_ele.PyPICobj.build_sparse_solver() if self.spacech_ele.flag_em_tracking: - print('Restoring PyPIC Ax, Ay and As state objects...') + print("Restoring PyPIC Ax, Ay and As state objects...") self.spacech_ele.state_Ax = self.spacech_ele.PyPICobj.get_state_object() self.spacech_ele.state_Ay = self.spacech_ele.PyPICobj.get_state_object() self.spacech_ele.state_As = self.spacech_ele.PyPICobj.get_state_object() - print('Done reload.') + print("Done reload.") return dict_state - def load_checkpoint(self, filename_simulation_checkpoint, load_from_folder='./'): - print(('Reloading from checkpoint: %s...' % (load_from_folder + filename_simulation_checkpoint))) - - i_checkp = int(filename_simulation_checkpoint.split('.pkl')[0].split('_')[-1]) - dict_state = self.load_state(filename_simulation_checkpoint, force_disable_save_simulation_state=False, filen_main_outp=None, load_from_folder=load_from_folder) + def load_checkpoint(self, filename_simulation_checkpoint, load_from_folder="./"): + print( + ( + "Reloading from checkpoint: %s..." + % (load_from_folder + filename_simulation_checkpoint) + ) + ) + + i_checkp = int(filename_simulation_checkpoint.split(".pkl")[0].split("_")[-1]) + dict_state = self.load_state( + filename_simulation_checkpoint, + force_disable_save_simulation_state=False, + filen_main_outp=None, + load_from_folder=load_from_folder, + ) for cloud in self.cloud_list: cloud.pyeclsaver.load_from_output(last_t=self.beamtim.tt_curr) cloud.pyeclsaver.i_checkp = i_checkp + 1 cloud.pyeclsaver.t_last_checkp = self.beamtim.tt_curr - cloud.pyeclsaver.t_last_En_hist = dict_state['t_last_En_hist'] + cloud.pyeclsaver.t_last_En_hist = dict_state["t_last_En_hist"] diff --git a/cloud_manager.py b/cloud_manager.py index 3908f88e..a6f1b9ba 100644 --- a/cloud_manager.py +++ b/cloud_manager.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Author and contact: Giovanni IADAROLA diff --git a/cross_ionization.py b/cross_ionization.py index 10349b25..040d0696 100644 --- a/cross_ionization.py +++ b/cross_ionization.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/default_input_parameters.py b/default_input_parameters.py index e2901cd9..6efbcd97 100644 --- a/default_input_parameters.py +++ b/default_input_parameters.py @@ -96,9 +96,10 @@ 'lifetime_hist_max': None, 'Dt_lifetime_hist':None, - 'Dh_electric_energy': None, + 'Dh_electric_energy': None, 'sparse_solver': 'scipy_slu', 'PyPICmode' : 'FiniteDifferences_ShortleyWeller', + 'flag_reinterp_fields_at_substeps': False, # Multigrid parameters 'f_telescope': None, diff --git a/doc/reference/reference.pdf b/doc/reference/reference.pdf index 37489e53..3c9cdd16 100644 Binary files a/doc/reference/reference.pdf and b/doc/reference/reference.pdf differ diff --git a/doc/reference/src/reference.tex b/doc/reference/src/reference.tex index 2be39f14..950297b9 100644 --- a/doc/reference/src/reference.tex +++ b/doc/reference/src/reference.tex @@ -205,7 +205,8 @@ \subsection{Simulation parameters} \textbf{sparse\_solver} & (optional) choices are 'klu', 'scipy\_slu' [default] \\\hline \textbf{PyPICmode} & (optional) Options are 'FiniteDifferences\_ShortleyWeller' [default], 'ShortleyWeller\_WithTelescopicGrids', 'FiniteDifferences\_Staircase' \\ \hline \textbf{Dh\_electric\_energy} & [m] (optional -- default = None) Grid size used to compute stored electric energy. If None is given, electric energy is not computed. \\ \hline - \textbf{flag\_em\_tracking} & (optional -- default = False) Forces generated by electrons are computed using electromagnetic potentials instead of quasi-electrostatic approximation (for more information see \url{https://indico.cern.ch/event/840676/contributions/3532754/}). + \textbf{flag\_em\_tracking} & (optional -- default = False) Forces generated by electrons are computed using electromagnetic potentials instead of quasi-electrostatic approximation (for more information see \url{https://indico.cern.ch/event/840676/contributions/3532754/}). \\ \hline + \textbf{flag\_reinterp\_fields\_at\_substeps} & (optional -- default = False) If true, field maps from beams and clouds are interpolated at each Boris substep. \\ \hline \end{longtable} diff --git a/dynamics_Boris_f2py.py b/dynamics_Boris_f2py.py index 6e0fdc93..ab803535 100644 --- a/dynamics_Boris_f2py.py +++ b/dynamics_Boris_f2py.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA @@ -131,6 +131,7 @@ def __init__(self, Dt, B0x, B0y, B0z, \ print("Tracker: Boris") + self.Dt = Dt self.N_sub_steps = N_sub_steps self.Dtt = Dt / float(N_sub_steps) self.B0x = B0x diff --git a/dynamics_Boris_multipole.py b/dynamics_Boris_multipole.py index 315b9078..e11444bb 100644 --- a/dynamics_Boris_multipole.py +++ b/dynamics_Boris_multipole.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/dynamics_dipole.py b/dynamics_dipole.py index 3722d91f..66a04a14 100644 --- a/dynamics_dipole.py +++ b/dynamics_dipole.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/dynamics_strong_B_generalized.py b/dynamics_strong_B_generalized.py index 482c6eaf..74eac7be 100644 --- a/dynamics_strong_B_generalized.py +++ b/dynamics_strong_B_generalized.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/electron_emission.py b/electron_emission.py index d03b9265..ec16c259 100644 --- a/electron_emission.py +++ b/electron_emission.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/gas_ionization_class.py b/gas_ionization_class.py index fe769a7d..512796c3 100644 --- a/gas_ionization_class.py +++ b/gas_ionization_class.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/gen_photoemission_class.py b/gen_photoemission_class.py index da47a951..ff23813e 100644 --- a/gen_photoemission_class.py +++ b/gen_photoemission_class.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/geom_impact_ellip.py b/geom_impact_ellip.py index 77f81c24..d9ef22f3 100644 --- a/geom_impact_ellip.py +++ b/geom_impact_ellip.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/geom_impact_poly.py b/geom_impact_poly.py index b7a67e9f..ba5bd9ab 100644 --- a/geom_impact_poly.py +++ b/geom_impact_poly.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/geom_impact_poly_fast_impact.py b/geom_impact_poly_fast_impact.py index 9cd6dbca..99699e50 100644 --- a/geom_impact_poly_fast_impact.py +++ b/geom_impact_poly_fast_impact.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/geom_impact_rect_fast_impact.py b/geom_impact_rect_fast_impact.py index e8efdeac..6c4ac99a 100644 --- a/geom_impact_rect_fast_impact.py +++ b/geom_impact_rect_fast_impact.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/h000_find_and_modify_preamble.py b/h000_find_and_modify_preamble.py index 69a4bcb4..462e4c64 100644 --- a/h000_find_and_modify_preamble.py +++ b/h000_find_and_modify_preamble.py @@ -16,7 +16,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/impact_management_class.py b/impact_management_class.py index ace7a18a..05faec91 100644 --- a/impact_management_class.py +++ b/impact_management_class.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/init.py b/init.py index 0637cccc..de373d12 100644 --- a/init.py +++ b/init.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA @@ -551,5 +551,6 @@ def read_input_files_and_init_components(pyecl_input_folder='./', skip_beam=Fals flag_multiple_clouds, cloud_list, cc.checkpoint_folder, - cross_ion + cross_ion, + cc.flag_reinterp_fields_at_substeps ) diff --git a/main.py b/main.py index 47b39582..6364c20b 100644 --- a/main.py +++ b/main.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/other/electron_oscillations/000_electron_oscillations.py b/other/electron_oscillations/000_electron_oscillations.py new file mode 100644 index 00000000..4b44b9e9 --- /dev/null +++ b/other/electron_oscillations/000_electron_oscillations.py @@ -0,0 +1,134 @@ +import time + +import numpy as np +from scipy.constants import e as qe +import matplotlib.pyplot as plt +from PyPARIS_sim_class import LHC_custom +from PyECLOUD import PyEC4PyHT +from PyHEADTAIL.particles.slicing import UniformBinSlicer + +import Simulation_parameters as pp + +plt.close('all') + +n_mp_slice = 25000 +intensity = 1.2e11 +epsn_x = 2.5e-6 +epsn_y = 2.5e-6 +sigma_z = 9.7e-2 +n_slices = 150 +n_mp_bunch = n_mp_slice*n_slices +long_unif = True + +################# +# Build machine # +################# + +machine = LHC_custom.LHC(n_segments=1, + machine_configuration='HLLHC-injection' + ) + +################ +# Make a bunch # +################ + +bunch = machine.generate_6D_Gaussian_bunch( + n_mp_bunch, intensity, epsn_x, epsn_y, sigma_z) +if long_unif: + bunch.z = 4*sigma_z*(np.random.rand(len(bunch.z))-0.5) + +########## +# Slicer # +########## + +slicer = UniformBinSlicer( + n_slices = n_slices, z_cuts=(-pp.z_cut, pp.z_cut)) + +######### +# Slice # +######### + +slices = bunch.extract_slices(slicer) +z_slices = np.array([ + ss.slice_info['z_bin_center'] for ss in slices]) +intensity_slices = np.array([ + ss.intensity for ss in slices]) + + +fig1 = plt.figure(1) +axl = fig1.add_subplot(111) +axl.plot(z_slices, intensity_slices) + +################## +# Make an ecloud # +################## + +inj_opt = machine.transverse_map.get_injection_optics() +beta_x_smooth = inj_opt['beta_x'] +beta_y_smooth = inj_opt['beta_y'] +sigma_x_smooth = np.sqrt(beta_x_smooth*epsn_x/machine.betagamma) +sigma_y_smooth = np.sqrt(beta_y_smooth*epsn_y/machine.betagamma) + +target_grid_arcs = { + 'x_min_target':-pp.target_size_internal_grid_sigma*sigma_x_smooth, + 'x_max_target':pp.target_size_internal_grid_sigma*sigma_x_smooth, + 'y_min_target':-pp.target_size_internal_grid_sigma*sigma_y_smooth, + 'y_max_target':pp.target_size_internal_grid_sigma*sigma_y_smooth, + 'Dh_target':pp.target_Dh_internal_grid_sigma*sigma_x_smooth} + +pp.N_mp_max_dip = 500000 +nel_mp_ref_0 = pp.init_unif_edens_dip*4*pp.x_aper*pp.y_aper/pp.N_MP_ele_init_dip + +ecloud = PyEC4PyHT.Ecloud(slice_by_slice_mode=True, + L_ecloud=1., slicer=None, + force_interp_at_substeps_interacting_slices=True, + Dt_ref=pp.Dt_ref, + pyecl_input_folder=pp.pyecl_input_folder, + chamb_type = pp.chamb_type, + x_aper=pp.x_aper, y_aper=pp.y_aper, + filename_chm=pp.filename_chm, + PyPICmode = pp.PyPICmode, + Dh_sc=pp.Dh_sc_ext, + N_min_Dh_main = pp.N_min_Dh_main, + f_telescope = pp.f_telescope, + N_nodes_discard = pp.N_nodes_discard, + target_grid = target_grid_arcs, + init_unif_edens_flag=pp.init_unif_edens_flag_dip, + init_unif_edens=pp.init_unif_edens_dip, + N_mp_max=pp.N_mp_max_dip, + nel_mp_ref_0=nel_mp_ref_0, + B_multip=[0.], + enable_kick_x = pp.enable_kick_x, + enable_kick_y = pp.enable_kick_y, + kick_mode_for_beam_field=False) + +sc = ecloud.beam_PyPIC_state.pic_internal +i_y0 = np.argmin(np.abs(sc.yg)) + + +mpe = ecloud.cloudsim.cloud_list[0].MP_e +mpe.x_mp[0] = .1e-3 +mpe.vx_mp[0] = 0. +mpe.y_mp[0] = 0. +mpe.vy_mp[0] = 0. +mpe.z_mp[0] = 0. +N_track = 1 +x_record = [] +field_record = [] + +t_start = time.mktime(time.localtime()) +for ss in slices[::-1]: + #ecloud.track(slices[int(n_slices/2)]) + ecloud.track(ss) + field_record.append(sc.efx[:, i_y0].copy()) + x_record.append(mpe.x_mp[:N_track].copy()) +x_record = np.array(x_record[::-1]) +t_end = time.mktime(time.localtime()) +print(f'Elapsed time {(t_end - t_start)}') + +fig2 = plt.figure(2) +axx = fig2.add_subplot(111) +axx.plot(z_slices, x_record, '.-') +axx.grid(True) +plt.show() + diff --git a/other/electron_oscillations/001_simple_euler.py b/other/electron_oscillations/001_simple_euler.py new file mode 100644 index 00000000..e59ffb85 --- /dev/null +++ b/other/electron_oscillations/001_simple_euler.py @@ -0,0 +1,32 @@ +import numpy as np + +N_steps = 100 +N_substeps = 10 +k_field = -0.05 + +Dt = 1. + +x =1e-3 +vx = 0. + +x_record = [] + +#for ii in range(N_steps): +# E = k_field*x +# vx += E*Dt +# x += vx*Dt +# x_record.append(x) + +for ii in range(N_steps): + E = k_field*x + for ssn in range(N_substeps): + vx += E*Dt/N_substeps + x += vx*Dt/N_substeps + x_record.append(x) + +import matplotlib.pyplot as plt +plt.close('all') +plt.plot(x_record) +plt.grid(True) +plt.show() + diff --git a/other/electron_oscillations/LHC_chm_ver.mat b/other/electron_oscillations/LHC_chm_ver.mat new file mode 100755 index 00000000..2a41ac7b Binary files /dev/null and b/other/electron_oscillations/LHC_chm_ver.mat differ diff --git a/other/electron_oscillations/Simulation_parameters.py b/other/electron_oscillations/Simulation_parameters.py new file mode 100644 index 00000000..d489a253 --- /dev/null +++ b/other/electron_oscillations/Simulation_parameters.py @@ -0,0 +1,149 @@ +from scipy.constants import c + +check_for_resubmit = False + +#################### +# Machine Settings # +#################### + +machine_configuration = 'LHC-collision' + +# # Use this part for optics from file +# # n_segments needs to be None if optics_pickle_file is specified +# optics_pickle_file = 'lhc2018_25cm_only_triplets_IR15_b1_optics.pkl' +# n_segments = None +# beta_x = None +# beta_y = None +# Q_x = None +# Q_y = None + +# # Use this part for smooth machine +optics_pickle_file = None +n_segments = 3 +beta_x = 400.0 +beta_y = 400.0 +Q_x = 62.27 +Q_y = 60.295 + +Qp_x = 0. +Qp_y = 0. + +octupole_knob = 0. + +V_RF = 12e6 + +n_non_parallelizable = 2 #rf and aperture + +# Transverse Damper Settings +enable_transverse_damper = False +dampingrate_x = 50. +dampingrate_y = 100. +if enable_transverse_damper: n_non_parallelizable += 1 + + +################### +# Beam Parameters # +################### + +bunch_from_file = None + +intensity = 1.2e+11 + +epsn_x = 2.5e-6 +epsn_y = 2.5e-6 + +sigma_z = 1.2e-9/4*c + +x_kick_in_sigmas = 0.1 +y_kick_in_sigmas = 0.1 + +# Numerical Parameters +n_slices = 100 +z_cut = 2.5e-9/2*c # For slicing +macroparticles_per_slice = 1000 +n_macroparticles = macroparticles_per_slice*n_slices + + +################# +# Stop Criteria # +################# + +# 1. Turns +N_turns = 10 # Per job +N_turns_target = 150 +# 2. Losses +sim_stop_frac = 0.9 +# 3. Emittance Growth +flag_check_emittance_growth = True +epsn_x_max_growth_fraction = 0.5 +epsn_y_max_growth_fraction = epsn_x_max_growth_fraction + + +###################### +# Footprint Settings # +###################### + +footprint_mode = False +n_macroparticles_for_footprint_map = 500000 +n_macroparticles_for_footprint_track = 5000 + + +#################### +# E-Cloud Settings # +#################### + +# General E-Cloud Settings +chamb_type = 'polyg' +x_aper = 2.300000e-02 +y_aper = 1.800000e-02 +filename_chm = 'LHC_chm_ver.mat' +Dt_ref = 5e-12 +pyecl_input_folder = './pyecloud_config' +sey = 1.30 + +# Transverse Multigrid Parameters +PyPICmode = 'ShortleyWeller_WithTelescopicGrids' +N_min_Dh_main = 10. +Dh_sc_ext = .8e-3 +f_telescope = 0.3 +N_nodes_discard = 5. +target_size_internal_grid_sigma = 10. +target_Dh_internal_grid_sigma = 0.2 +custom_target_grid_arcs = None + +# # Uncomment for custom grid +# custom_target_grid_arcs = { +# 'x_min_target': -3e-3, +# 'x_max_target': 3e-3, +# 'y_min_target': -3.1e-3, +# 'y_max_target': 3.1e-3, +# 'Dh_target': 7e-5} + +# Enable Kicks Different Planes +enable_kick_x = True +enable_kick_y = True + +# Dedicated Dipole E-Cloud Settings +enable_arc_dip = True +fraction_device_dip = 0.65 +init_unif_edens_flag_dip = 1 +init_unif_edens_dip = 1.000000e+12 +N_MP_ele_init_dip = 50000 +N_mp_max_dip = N_MP_ele_init_dip*4 +B_multip_dip = [8.33] #T + +# Dedicated Quadrupole E-Cloud Settings +enable_arc_quad = False +fraction_device_quad = 26.000000e-02 #7.000000e-02 +N_mp_max_quad = 2000000 +B_multip_quad = [0., 188.2] #T +folder_path = '../../LHC_ecloud_distrib_quads/' +filename_state = 'combined_distribution_sey%.2f_%.1fe11ppb_7tev.mat'%(sey,intensity/1e11) +filename_init_MP_state_quad = folder_path + filename_state + +# Dedicated Kick Element Settings +enable_eclouds_at_kick_elements = False +path_buildup_simulations_kick_elements = '/home/kparasch/workspace/Triplets/ec_headtail_triplets/simulations_PyECLOUD/!!!NAME!!!_sey1.35' +name_MP_state_file_kick_elements = 'MP_state_9.mat' +orbit_factor = 6.250000e-01 + diff --git a/other/electron_oscillations/pyecloud_config/machine_parameters.input b/other/electron_oscillations/pyecloud_config/machine_parameters.input new file mode 100644 index 00000000..0054aec7 --- /dev/null +++ b/other/electron_oscillations/pyecloud_config/machine_parameters.input @@ -0,0 +1,11 @@ +# MACHINE PARAMETERS + +chamb_type = None #redefine in the script + +photoelectron_angle_distribution = 'cosine_3D' + +track_method= 'BorisMultipole' +N_sub_steps = 5 +B_multip = [0.] + + diff --git a/other/electron_oscillations/pyecloud_config/secondary_emission_parameters.input b/other/electron_oscillations/pyecloud_config/secondary_emission_parameters.input new file mode 100644 index 00000000..ef0930c6 --- /dev/null +++ b/other/electron_oscillations/pyecloud_config/secondary_emission_parameters.input @@ -0,0 +1,17 @@ +# secondary emission model +Emax=332.; +del_max = 0. +R0 = 0. +switch_model=0 + +# hilleret model for sec. emission +E_th=35.; +sigmafit =1.0828; +mufit = 1.6636; + +switch_no_increase_energy=0 +thresh_low_energy=-1 + +scrub_en_th=20.#eV + +secondary_angle_distribution = 'cosine_3D' diff --git a/other/electron_oscillations/pyecloud_config/simulation_parameters.input b/other/electron_oscillations/pyecloud_config/simulation_parameters.input new file mode 100644 index 00000000..11b82691 --- /dev/null +++ b/other/electron_oscillations/pyecloud_config/simulation_parameters.input @@ -0,0 +1,63 @@ +# SIMULATION PARAMETERS + +machine_param_file='machine_parameters.input' +secondary_emission_parameters_file='secondary_emission_parameters.input' +beam_parameters_file='beam.beam' + +logfile_path = 'logfile.txt' +progress_path = 'progress' +stopfile = 'stop' + +Dt = None +t_end = None #s (no effect if log. profile is imported from file) + +#import numpy as np +#dec_fact_out = int(np.round(5 * 25e-12/Dt)) + +lam_th = None #e-/m +Dx_hist = 1e-3 #m +r_center = 1e-3 #m + + +Dt_En_hist = 25e-9 #s +Nbin_En_hist= 2000 +En_hist_max= 2e3 #eV + +t_ion=100.; #s + +N_mp_max=250000; #size of allocated vectors + +#Regen parameters + +N_mp_regen=250000; +N_mp_regen_low=5000; +N_mp_after_regen=10000; +t_ON_regen_low=10. +fact_split=1.5; +fact_clean=1e-6; +regen_hist_cut = 1.e-4 + +N_mp_soft_regen = 75000 +N_mp_after_soft_regen = 25000 + +nel_mp_ref_0 = None #redefine in the script + +# Number of bins +Nx_regen=51;#it must be odd! +Ny_regen=51;#it must be odd! +Nvx_regen=51;#it must be odd! +Nvy_regen=101;#it must be odd! +Nvz_regen=51;#it must be odd! + + +#Sp_ch params +Dt_sc = None +Dh_sc = None +t_sc_ON=0e-9; #s +#PyPICmode = 'FFT_OpenBoundary' +sparse_solver = 'klu' + +flag_movie = 0 #1/0 +flag_sc_movie = 0 #1/0 + +save_mp_state_time_file = -1 diff --git a/proc_video3.py b/proc_video3.py index fb294cfb..3e8f1a38 100644 --- a/proc_video3.py +++ b/proc_video3.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/pyecloud_saver.py b/pyecloud_saver.py index 9ed09f1c..95ff6e59 100644 --- a/pyecloud_saver.py +++ b/pyecloud_saver.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA @@ -96,7 +96,7 @@ def __init__(self, logfile_path): if self.logfile_path is not None: with open(self.logfile_path, 'w') as flog: - flog.write('PyECLOUD Version 8.2.0\n') + flog.write('PyECLOUD Version 8.3.0\n') flog.write('%s\n' % git_hash) flog.write('%s\n' % git_branch) flog.write('Simulation started on %s\n' % timestr) diff --git a/sec_emission_model_ECLOUD.py b/sec_emission_model_ECLOUD.py index 7717683e..e7969c5e 100644 --- a/sec_emission_model_ECLOUD.py +++ b/sec_emission_model_ECLOUD.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_ECLOUD_nunif.py b/sec_emission_model_ECLOUD_nunif.py index 56409cda..a5975b9d 100644 --- a/sec_emission_model_ECLOUD_nunif.py +++ b/sec_emission_model_ECLOUD_nunif.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_accurate_low_ene.py b/sec_emission_model_accurate_low_ene.py index 0c05db1a..13df3917 100644 --- a/sec_emission_model_accurate_low_ene.py +++ b/sec_emission_model_accurate_low_ene.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_cos_low_ener.py b/sec_emission_model_cos_low_ener.py index 6e8d44a2..12f299f4 100644 --- a/sec_emission_model_cos_low_ener.py +++ b/sec_emission_model_cos_low_ener.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_flat_low_ener.py b/sec_emission_model_flat_low_ener.py index 88f73b23..ace40539 100644 --- a/sec_emission_model_flat_low_ener.py +++ b/sec_emission_model_flat_low_ener.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_from_file.py b/sec_emission_model_from_file.py index ead56ef9..fa2b9a51 100644 --- a/sec_emission_model_from_file.py +++ b/sec_emission_model_from_file.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_furman_pivi.py b/sec_emission_model_furman_pivi.py index c09fa1a6..21c0deb7 100644 --- a/sec_emission_model_furman_pivi.py +++ b/sec_emission_model_furman_pivi.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sec_emission_model_perfect_absorber.py b/sec_emission_model_perfect_absorber.py index 13c89ed8..4eda3371 100644 --- a/sec_emission_model_perfect_absorber.py +++ b/sec_emission_model_perfect_absorber.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/sincc_cosincc_cubsincc.py b/sincc_cosincc_cubsincc.py index 1d0a09cc..4fc83064 100644 --- a/sincc_cosincc_cubsincc.py +++ b/sincc_cosincc_cubsincc.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/space_charge_class.py b/space_charge_class.py index 02815f97..a2199c4d 100644 --- a/space_charge_class.py +++ b/space_charge_class.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/space_charge_class_electromagnetic.py b/space_charge_class_electromagnetic.py index 2a01e2bd..890b7c07 100644 --- a/space_charge_class_electromagnetic.py +++ b/space_charge_class_electromagnetic.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyECLOUD Version 8.2.0 +# PyECLOUD Version 8.3.0 # # # Main author: Giovanni IADAROLA diff --git a/testing/tests_buildup/002_run_all_tests.py b/testing/tests_buildup/002_run_all_tests.py index 8f8e8ed6..b11675ed 100644 --- a/testing/tests_buildup/002_run_all_tests.py +++ b/testing/tests_buildup/002_run_all_tests.py @@ -34,7 +34,8 @@ 'LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_nonuniftime', 'LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_checkpoint', 'Rectangular_Dip_450GeV_sey1.60_1.1e11ppb_furman_pivi', - 'LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_em_tracking' + 'LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_em_tracking', + 'LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps', ] if args.all: diff --git a/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/LHC_chm_ver.mat b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/LHC_chm_ver.mat new file mode 100755 index 00000000..2a41ac7b Binary files /dev/null and b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/LHC_chm_ver.mat differ diff --git a/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/Pyecltest_angle3D_ref.mat b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/Pyecltest_angle3D_ref.mat new file mode 100644 index 00000000..dd37e206 Binary files /dev/null and b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/Pyecltest_angle3D_ref.mat differ diff --git a/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/beam.beam b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/beam.beam new file mode 100755 index 00000000..f24483bd --- /dev/null +++ b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/beam.beam @@ -0,0 +1,28 @@ +# BEAM PARAMETERS +energy_eV = 4.500000e+11 +nemittx = 2.500000e-06 +nemitty = nemittx #m +Dp_p = 0.000000e+00 + + +beam_field_file = 'computeFD' +Dh_beam_field = 0.0001 +Nx = None +Ny = None +nimag = None + + +b_spac = 25e-9 #s (to be specified also if you load the profile from file - + # it is used as period for clean and save) +fact_beam = 2.500000e+11 +coast_dens = 0. #protons per meter + +flag_bunched_beam = 1 # 1: bunched beam 0:load profile from file + +# to be filled in case of bunched beam +sigmaz = 1.000000e-09/4.*299792458. +t_offs = 2.5e-9 +filling_pattern_file = 1*(30*[1.]+5*[0]) + +# to be filled in case of longitudinal profile from file +beam_long_prof_file = -1 diff --git a/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/machine_parameters.input b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/machine_parameters.input new file mode 100755 index 00000000..558c4750 --- /dev/null +++ b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/machine_parameters.input @@ -0,0 +1,37 @@ +# MACHINE PARAMETERS + +chamb_type = 'polyg_cython' + +x_aper = 2.300000e-02 +y_aper = 1.800000e-02 + +filename_chm = 'LHC_chm_ver.mat' + + + + +# Choose track_method= 'StrongBdip' / 'StrongBgen'/ 'Boris' +track_method= 'BorisMultipole' +N_sub_steps = 5 +B_multip = [0.53549999999999998] + + + + +betafx = 85.00 +betafy = 90.00 + + +# uniform initial distrib +init_unif_flag=1 +Nel_init_unif=2e8 +E_init_unif=0.01; #in eV +x_max_init_unif =x_aper +x_min_init_unif =-x_aper +y_max_init_unif =y_aper +y_min_init_unif =-y_aper + + + + +photoelectron_angle_distribution = 'cosine_3D' diff --git a/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/secondary_emission_parameters.input b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/secondary_emission_parameters.input new file mode 100755 index 00000000..7fbaeb23 --- /dev/null +++ b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/secondary_emission_parameters.input @@ -0,0 +1,18 @@ +# secondary emission model +Emax=332.; +del_max = 1.700000 +R0 = 0.7 +switch_model=0 + +# hilleret model for sec. emission +E_th=35.; +sigmafit =1.0828; +mufit = 1.6636; + +switch_no_increase_energy=0 +thresh_low_energy=-1 + +scrub_en_th=20.#eV + + +secondary_angle_distribution = 'cosine_3D' diff --git a/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/simulation_parameters.input b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/simulation_parameters.input new file mode 100755 index 00000000..f5fd8c12 --- /dev/null +++ b/testing/tests_buildup/LHC_ArcDipReal_450GeV_sey1.70_2.5e11ppb_bl_1.00ns_reinterp_at_substeps/simulation_parameters.input @@ -0,0 +1,64 @@ +# SIMULATION PARAMETERS + +machine_param_file='machine_parameters.input' +secondary_emission_parameters_file='secondary_emission_parameters.input' +beam_parameters_file='beam.beam' + +logfile_path = './logfile.txt' +progress_path = './progress.txt' +stopfile = 'stop' + +Dt = 2.500000e-11 +t_end=1e-9; #s (no effect if log. profile is imported from file) + +import numpy as np +dec_fact_out = int(np.round(5 * 25e-12/Dt)) + +lam_th=1.e2 #e-/m +Dx_hist=1.e-3 #m +r_center=1.e-3 #m + + +Dt_En_hist = 25e-9 #s +Nbin_En_hist= 250 +En_hist_max= 5000. #eV + +t_ion=100.; #s + +N_mp_max=250000; #size of allocated vectors + +#Regen parameters + +N_mp_regen=200000; +N_mp_regen_low=5000; +N_mp_after_regen=10000; +t_ON_regen_low=10. +fact_split=1.5; +fact_clean=1e-6; +regen_hist_cut = 1.e-4 + +N_mp_soft_regen = 60000 +N_mp_after_soft_regen = 20000 + +nel_mp_ref_0= 2e8/(0.7*N_mp_soft_regen) #e-/m + + +# Number of bins +Nx_regen=51;#it must be odd! +Ny_regen=51;#it must be odd! +Nvx_regen=51;#it must be odd! +Nvy_regen=101;#it must be odd! +Nvz_regen=51;#it must be odd! + + +#Sp_ch params +Dt_sc = .5e-9 +Dh_sc = .3e-3 +t_sc_ON=0e-9; #s +sparse_solver = 'klu' +flag_reinterp_fields_at_substeps = True + +flag_movie = 0 #1/0 +flag_sc_movie = 0 #1/0 + +save_mp_state_time_file = -1