diff --git a/Bassetti_Erskine.py b/Bassetti_Erskine.py index cf8b6d3a..97a286de 100644 --- a/Bassetti_Erskine.py +++ b/Bassetti_Erskine.py @@ -8,7 +8,7 @@ # This file is part of the code # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/FFT_OpenBoundary.py b/FFT_OpenBoundary.py index e0361f44..e0cdc6a0 100644 --- a/FFT_OpenBoundary.py +++ b/FFT_OpenBoundary.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/FFT_OpenBoundary_SquareGrid.py b/FFT_OpenBoundary_SquareGrid.py index 48e16abb..30e70bc3 100644 --- a/FFT_OpenBoundary_SquareGrid.py +++ b/FFT_OpenBoundary_SquareGrid.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/FFT_PEC_Boundary_SquareGrid.py b/FFT_PEC_Boundary_SquareGrid.py index c1e7ab74..446e9954 100644 --- a/FFT_PEC_Boundary_SquareGrid.py +++ b/FFT_PEC_Boundary_SquareGrid.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/FiniteDifferences_ShortleyWeller_SquareGrid.py b/FiniteDifferences_ShortleyWeller_SquareGrid.py index 9d94cee0..20646aa1 100644 --- a/FiniteDifferences_ShortleyWeller_SquareGrid.py +++ b/FiniteDifferences_ShortleyWeller_SquareGrid.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/FiniteDifferences_ShortleyWeller_SquareGrid_extrapolation.py b/FiniteDifferences_ShortleyWeller_SquareGrid_extrapolation.py index 48a5575b..4737ab08 100644 --- a/FiniteDifferences_ShortleyWeller_SquareGrid_extrapolation.py +++ b/FiniteDifferences_ShortleyWeller_SquareGrid_extrapolation.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/FiniteDifferences_Staircase_SquareGrid.py b/FiniteDifferences_Staircase_SquareGrid.py index 9f343b7c..7e7f22cb 100644 --- a/FiniteDifferences_Staircase_SquareGrid.py +++ b/FiniteDifferences_Staircase_SquareGrid.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/PyPIC_Scatter_Gather.py b/PyPIC_Scatter_Gather.py index 3c9f4b19..831e5ccd 100644 --- a/PyPIC_Scatter_Gather.py +++ b/PyPIC_Scatter_Gather.py @@ -8,7 +8,7 @@ # This file is part of the code: # # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA @@ -66,7 +66,7 @@ class PyPIC_Scatter_Gather(object): def __init__(self, x_aper=None, y_aper=None, dx=None, dy=None, xg=None, yg=None, x_min=None, x_max=None, y_min=None, y_max=None, *args, **kwargs): - print('PyPIC Version 2.4.4') + print('PyPIC Version 2.4.5') if xg is not None and yg is not None: assert(x_aper is None and y_aper is None and dx is None and dy is None) diff --git a/change_version_number.py b/change_version_number.py index 6c27779d..c1439231 100644 --- a/change_version_number.py +++ b/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('PyPIC Version 2.4.4', 'PyPIC Version 2.4.4') + content=content.replace('PyPIC Version 2.4.5', 'PyPIC Version 2.4.5') with open(filename,'w') as fid: fid.write(content) diff --git a/geom_impact_poly.py b/geom_impact_poly.py index f562a55d..21c3ab81 100644 --- a/geom_impact_poly.py +++ b/geom_impact_poly.py @@ -7,7 +7,7 @@ # # This file is part of the code: # -# PyPIC Version 2.4.4 +# PyPIC Version 2.4.5 # # # Author and contact: Giovanni IADAROLA diff --git a/temp_test_states/000_test_states_FDSW.py b/temp_test_states/000_test_states_FDSW.py deleted file mode 100644 index b102373a..00000000 --- a/temp_test_states/000_test_states_FDSW.py +++ /dev/null @@ -1,124 +0,0 @@ -import sys, os -BIN=os.path.expanduser('../../') -sys.path.append(BIN) -BIN=os.path.expanduser('../../PyHEADTAIL/testing/script-tests/') -sys.path.append(BIN) -from LHC import LHC -import PyPIC.geom_impact_ellip as ell -import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW -import numpy as np -import pylab as pl -import mystyle as ms -from scipy.constants import e, epsilon_0 - - -qe = e -eps0 = epsilon_0 - -# chamber parameters -x_aper = 25e-3 -y_aper = 25e-3 -Dh_single = 0.5e-3 -#~ # machine parameters -optics_mode = 'smooth' -n_segments=1 -# beam parameters -n_macroparticles=1000000 - - -#LHC -machine_configuration='6.5_TeV_collision_tunes' -intensity=1.2e11 -epsn_x=.5e-6 -epsn_y=3e-6 -sigma_z=7e-2 -machine = LHC(machine_configuration = machine_configuration, optics_mode = optics_mode, n_segments = n_segments) - - -# build chamber -chamber = ell.ellip_cham_geom_object(x_aper = x_aper, y_aper = y_aper) -Vx, Vy = chamber.points_on_boundary(N_points=200) - - -# generate beam for singlegrid -bunch = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state1 -bunch1 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 2*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state2 -bunch2 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 3*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state3 -bunch3 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 5*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -#scatter and solve singlegrid -pic_singlegrid = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_single) -pic_singlegrid.scatter(bunch.x, bunch.y, bunch.particlenumber_per_mp+bunch.y*0., charge=qe) -pic_singlegrid.solve() - -#states -state1 = pic_singlegrid.get_state_object() -state2 = pic_singlegrid.get_state_object() -state3 = pic_singlegrid.get_state_object() - -#scatter states -state1.scatter(bunch1.x, bunch1.y, bunch1.particlenumber_per_mp+bunch1.y*0., charge=qe) -state2.scatter(bunch2.x, bunch2.y, bunch2.particlenumber_per_mp+bunch2.y*0., charge=qe) -state3.scatter(bunch3.x, bunch3.y, bunch3.particlenumber_per_mp+bunch3.y*0., charge=qe) - -#solve states -pic_singlegrid.solve_states([state1, state2, state3]) - - - - -#plot electric field for each state and singlegrid -pl.close('all') - -# prepare probes -theta_probes=np.linspace(0., 2*np.pi, 100) -r_probes= 0.2e-3 -x_probes = r_probes*np.cos(theta_probes) -y_probes = r_probes*np.sin(theta_probes) - -# get field at probes -Ex_singlegrid, Ey_singlegrid = pic_singlegrid.gather(x_probes, y_probes) -Ex_state1, Ey_state1 = state1.gather(x_probes, y_probes) -Ex_state2, Ey_state2 = state2.gather(x_probes, y_probes) -Ex_state3, Ey_state3 = state3.gather(x_probes, y_probes) - - - -#~ #plot at probes -pl.close('all') -ms.mystyle_arial(fontsz=12) -pl.figure(4, figsize=(8,6)).patch.set_facecolor('w') -sp1=pl.subplot(2,1,1) -pl.plot(theta_probes, Ex_singlegrid, '--k', label = 'Singlegrid (I)') -pl.plot(theta_probes, Ex_state1, '.-m', label = 'State 1 (2I)') -pl.plot(theta_probes, Ex_state2, '.-y', label = 'State 2 (3I)') -pl.plot(theta_probes, Ex_state3, '.-g', label = 'State 3 (5I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ex [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.grid('on') -pl.subplot(2,1,2, sharex=sp1) -pl.plot(theta_probes, Ey_singlegrid, '--k', label = 'Singlegrid (I)') -pl.plot(theta_probes, Ey_state1, '.-m', label = 'State 1 (2I)') -pl.plot(theta_probes, Ey_state2, '.-y', label = 'State 2 (3I)') -pl.plot(theta_probes, Ey_state3, '.-g', label = 'State 3 (5I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ey [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.suptitle(' Test states for FiniteDifferences_ShortleyWeller_SquareGrid ') -pl.grid('on') -pl.show() diff --git a/temp_test_states/001_test_states_FDSC.py b/temp_test_states/001_test_states_FDSC.py deleted file mode 100644 index 6b5bfc40..00000000 --- a/temp_test_states/001_test_states_FDSC.py +++ /dev/null @@ -1,121 +0,0 @@ -import sys, os -BIN=os.path.expanduser('../../') -sys.path.append(BIN) -BIN=os.path.expanduser('../../PyHEADTAIL/testing/script-tests/') -sys.path.append(BIN) -from LHC import LHC -import PyPIC.geom_impact_ellip as ell -import PyPIC.FiniteDifferences_Staircase_SquareGrid as PIC_FDSC -import numpy as np -import pylab as pl -import mystyle as ms -from scipy.constants import e, epsilon_0 - - -qe = e -eps0 = epsilon_0 - -# chamber parameters -x_aper = 25e-3 -y_aper = 25e-3 -Dh_single = 0.5e-3 -#~ # machine parameters -optics_mode = 'smooth' -n_segments=1 -# beam parameters -n_macroparticles=1000000 - - -#LHC -machine_configuration='6.5_TeV_collision_tunes' -intensity=1.2e11 -epsn_x=.5e-6 -epsn_y=3e-6 -sigma_z=7e-2 -machine = LHC(machine_configuration = machine_configuration, optics_mode = optics_mode, n_segments = n_segments) - - -# build chamber -chamber = ell.ellip_cham_geom_object(x_aper = x_aper, y_aper = y_aper) -Vx, Vy = chamber.points_on_boundary(N_points=200) - - -# generate beam for singlegrid -bunch = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state1 -bunch1 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 2*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state2 -bunch2 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 3*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state3 -bunch3 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 5*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -#scatter and solve singlegrid -pic_singlegrid = PIC_FDSC.FiniteDifferences_Staircase_SquareGrid(chamb = chamber, Dh = Dh_single) -pic_singlegrid.scatter(bunch.x, bunch.y, bunch.particlenumber_per_mp+bunch.y*0., charge=qe) -pic_singlegrid.solve() - -#states -state1 = pic_singlegrid.get_state_object() -state2 = pic_singlegrid.get_state_object() -state3 = pic_singlegrid.get_state_object() - -#scatter states -state1.scatter(bunch1.x, bunch1.y, bunch1.particlenumber_per_mp+bunch1.y*0., charge=qe) -state2.scatter(bunch2.x, bunch2.y, bunch2.particlenumber_per_mp+bunch2.y*0., charge=qe) -state3.scatter(bunch3.x, bunch3.y, bunch3.particlenumber_per_mp+bunch3.y*0., charge=qe) - -#solve states -pic_singlegrid.solve_states([state1, state2, state3]) - - -#plot electric field for each state and singlegrid -pl.close('all') - -#prepare probes -theta_probes=np.linspace(0., 2*np.pi, 100) -r_probes= 0.2e-3 -x_probes = r_probes*np.cos(theta_probes) -y_probes = r_probes*np.sin(theta_probes) - -# get field at probes -Ex_singlegrid, Ey_singlegrid = pic_singlegrid.gather(x_probes, y_probes) -Ex_state1, Ey_state1 = state1.gather(x_probes, y_probes) -Ex_state2, Ey_state2 = state2.gather(x_probes, y_probes) -Ex_state3, Ey_state3 = state3.gather(x_probes, y_probes) - - -#plot at probes -pl.close('all') -ms.mystyle_arial(fontsz=12) -pl.figure(4, figsize=(8,6)).patch.set_facecolor('w') -sp1=pl.subplot(2,1,1) -pl.plot(theta_probes, Ex_singlegrid, '--k', label = 'Singlegrid (I)') -pl.plot(theta_probes, Ex_state1, '.-m', label = 'State 1 (2I)') -pl.plot(theta_probes, Ex_state2, '.-y', label = 'State 2 (3I)') -pl.plot(theta_probes, Ex_state3, '.-g', label = 'State 3 (5I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ex [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.grid('on') -pl.subplot(2,1,2, sharex=sp1) -pl.plot(theta_probes, Ey_singlegrid, '--k', label = 'Singlegrid (I)') -pl.plot(theta_probes, Ey_state1, '.-m', label = 'State 1 (2I)') -pl.plot(theta_probes, Ey_state2, '.-y', label = 'State 2 (3I)') -pl.plot(theta_probes, Ey_state3, '.-g', label = 'State 3 (5I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ey [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.suptitle('Test states for FiniteDifferences_Staircase_SquareGrid') -pl.grid('on') -pl.show() diff --git a/temp_test_states/002_test_states_InternalGrid.py b/temp_test_states/002_test_states_InternalGrid.py deleted file mode 100644 index a67b5d1d..00000000 --- a/temp_test_states/002_test_states_InternalGrid.py +++ /dev/null @@ -1,104 +0,0 @@ -import sys, os -BIN=os.path.expanduser('../../') -sys.path.append(BIN) -BIN=os.path.expanduser('../../PyHEADTAIL/testing/script-tests/') -sys.path.append(BIN) -from LHC import LHC -from PyPIC.MultiGrid import AddInternalGrid -import PyPIC.geom_impact_ellip as ell -import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW -import numpy as np -import pylab as pl -import mystyle as ms -from scipy.constants import e, epsilon_0 - - -qe = e -eps0 = epsilon_0 - -# chamber parameters -x_aper = 25e-3 -y_aper = 25e-3 -Dh_single = 0.5e-3 -#~ # machine parameters -optics_mode = 'smooth' -n_segments=1 -# beam parameters -n_macroparticles=1000000 - - -#LHC -machine_configuration='6.5_TeV_collision_tunes' -intensity=1.2e11 -epsn_x=.5e-6 -epsn_y=3e-6 -sigma_z=7e-2 -machine = LHC(machine_configuration = machine_configuration, optics_mode = optics_mode, n_segments = n_segments) - - -# build chamber -chamber = ell.ellip_cham_geom_object(x_aper = x_aper, y_aper = y_aper) -Vx, Vy = chamber.points_on_boundary(N_points=200) - -# generate beam for dualgrid -bunch_dual = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam for state -bunch_state = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 2*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -#create, scatter and solve dualgrid -pic_singlegrid = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_single) - -pic_dualgrid = AddInternalGrid(pic_singlegrid, x_min_internal = -2e-3, x_max_internal = 2e-3, y_min_internal = -1e-3, y_max_internal = 1e-3, Dh_internal = 0.2e-4, N_nodes_discard = 3, - sparse_solver = 'PyKLU', include_solver = True) - -pic_dualgrid.scatter(bunch_dual.x, bunch_dual.y, bunch_dual.particlenumber_per_mp+bunch_dual.y*0., charge=qe) -pic_dualgrid.solve() - -#state -state1 = pic_dualgrid.get_state_object() -#scatter state -state1.scatter(bunch_state.x, bunch_state.y, bunch_state.particlenumber_per_mp+bunch_state.y*0., charge=qe) -#solve state -pic_dualgrid.solve_states([state1]) - - -#plot electric field for each state and singlegrid -pl.close('all') -#~ #prepare probes -theta_probes=np.linspace(0., 2*np.pi, 100) -r_probes= 0.2e-3 -x_probes = r_probes*np.cos(theta_probes) -y_probes = r_probes*np.sin(theta_probes) - -# get field at probes -Ex_dualgrid, Ey_dualgrid = pic_dualgrid.gather(x_probes, y_probes) -Ex_state1, Ey_state1 = state1.gather(x_probes, y_probes) - - -#plot at probes -pl.close('all') -ms.mystyle_arial(fontsz=12) -pl.figure(4, figsize=(8,6)).patch.set_facecolor('w') -sp1=pl.subplot(2,1,1) -pl.plot(theta_probes, Ex_dualgrid, '--k', label = 'Dualgrid (I)') -pl.plot(theta_probes, Ex_state1, '.-m', label = 'State 1 (2I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ex [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.grid('on') -pl.subplot(2,1,2, sharex=sp1) -pl.plot(theta_probes, Ey_dualgrid, '--k', label = 'Dualgrid (I)') -pl.plot(theta_probes, Ey_state1, '.-m', label = 'State 1 (2I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ey [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.suptitle('Test states for AddInternalGrid') -pl.grid('on') -pl.show() diff --git a/temp_test_states/003_test_states_MultiGrid.py b/temp_test_states/003_test_states_MultiGrid.py deleted file mode 100644 index b9af2510..00000000 --- a/temp_test_states/003_test_states_MultiGrid.py +++ /dev/null @@ -1,136 +0,0 @@ -import sys, os -BIN=os.path.expanduser('../../') -sys.path.append(BIN) -BIN=os.path.expanduser('../../PyHEADTAIL/testing/script-tests/') -sys.path.append(BIN) -from LHC import LHC -import PyPIC.geom_impact_ellip as ell -from PyPIC.MultiGrid import AddTelescopicGrids -import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW -import numpy as np -import pylab as pl -import mystyle as ms -from scipy.constants import e, epsilon_0 - - -qe = e -eps0 = epsilon_0 - -# chamber parameters -x_aper = 25e-3 -y_aper = 25e-3 -Dh_single = 0.5e-3 -#~ # machine parameters -optics_mode = 'smooth' -n_segments=1 -# beam parameters -n_macroparticles=1000000 - - -#LHC -machine_configuration='6.5_TeV_collision_tunes' -intensity=1.2e11 -epsn_x=.5e-6 -epsn_y=3e-6 -sigma_z=7e-2 -machine = LHC(machine_configuration = machine_configuration, optics_mode = optics_mode, n_segments = n_segments) - - -# build chamber -chamber = ell.ellip_cham_geom_object(x_aper = x_aper, y_aper = y_aper) -Vx, Vy = chamber.points_on_boundary(N_points=200) - - -# generate beam for multigrid -bunch = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - - -# generate beam 1 for state1 -bunch1 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 2*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - - -# generate beam 2 for state 2 -bunch2 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 3*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - -# generate beam 3 for state3 -bunch3 = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = 5*intensity, - epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) - - -# Multi grid parameters -Sx_target = 10*bunch.sigma_x() -Sy_target = 10*bunch.sigma_y() -Dh_target = 0.5*bunch.sigma_x() -sparse_solver = 'PyKLU' - -pic_singlegrid = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_single) - -# build telescope -pic_multigrid = AddTelescopicGrids(pic_main = pic_singlegrid, f_telescope = 0.3, - target_grid = {'x_min_target':-Sx_target/2., 'x_max_target':Sx_target/2.,'y_min_target':-Sy_target/2.,'y_max_target':Sy_target/2.,'Dh_target':Dh_target}, - N_nodes_discard = 3., N_min_Dh_main = 10, sparse_solver=sparse_solver) -#scatter and solve multigrid -pic_multigrid.scatter(bunch.x, bunch.y, bunch.particlenumber_per_mp+bunch.y*0., charge=qe) -pic_multigrid.solve() - -#states -state1 = pic_multigrid.get_state_object() -state2 = pic_multigrid.get_state_object() -state3 = pic_multigrid.get_state_object() - -#scatter states -state1.scatter(bunch1.x, bunch1.y, bunch1.particlenumber_per_mp+bunch1.y*0., charge=qe) -state2.scatter(bunch2.x, bunch2.y, bunch2.particlenumber_per_mp+bunch2.y*0., charge=qe) -state3.scatter(bunch3.x, bunch3.y, bunch3.particlenumber_per_mp+bunch3.y*0., charge=qe) -#solve states -pic_multigrid.solve_states([state1, state2, state3]) - - - - -#plot electric field for each state and singlegrid -pl.close('all') - -# prepare probes -theta_probes=np.linspace(0., 2*np.pi, 100) -r_probes= 0.2e-3 -x_probes = r_probes*np.cos(theta_probes) -y_probes = r_probes*np.sin(theta_probes) - -# get field at probes -Ex_multigrid, Ey_multigrid = pic_multigrid.gather(x_probes, y_probes) -Ex_state1, Ey_state1 = state1.gather(x_probes, y_probes) -Ex_state2, Ey_state2 = state2.gather(x_probes, y_probes) -Ex_state3, Ey_state3 = state3.gather(x_probes, y_probes) - -#~ #plot at probes -pl.close('all') -ms.mystyle_arial(fontsz=12) -pl.figure(4, figsize=(8,6)).patch.set_facecolor('w') -sp1=pl.subplot(2,1,1) -pl.plot(theta_probes, Ex_multigrid, '--k', label = 'Multigrid (I)') -pl.plot(theta_probes, Ex_state1, '.-m', label = 'State 1 (2I)') -pl.plot(theta_probes, Ex_state2, '.-y', label = 'State 2 (3I)') -pl.plot(theta_probes, Ex_state3, '.-g', label = 'State 3 (5I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ex [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.grid('on') -pl.subplot(2,1,2, sharex=sp1) -pl.plot(theta_probes, Ey_multigrid, '--k', label = 'Multigrid (I)') -pl.plot(theta_probes, Ey_state1, '.-m', label = 'State 1 (2I)') -pl.plot(theta_probes, Ey_state2, '.-y', label = 'State 2 (3I)') -pl.plot(theta_probes, Ey_state3, '.-g', label = 'State 3 (5I)') -pl.xlabel('theta[deg]') -pl.ylabel('Ey [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') -pl.legend(loc = 'best') -pl.suptitle('Test states for MultiGrid') -pl.grid('on') -pl.show() diff --git a/temp_test_states/mystyle.py b/temp_test_states/mystyle.py deleted file mode 100644 index 73e2d51d..00000000 --- a/temp_test_states/mystyle.py +++ /dev/null @@ -1,27 +0,0 @@ -from matplotlib import rc, rcdefaults -import pylab as pl - - - -def mystyle(fontsz=10): - - font = {#'family' : 'normal', - #'weight' : 'bold', - 'size' : fontsz} - print(fontsz) - rcdefaults() - rc('font', **font) - -def mystyle_arial(fontsz=10, dist_tick_lab=7): - - rcdefaults() - rc('font',**{'family':'sans-serif','sans-serif':['arial'], 'size':fontsz}) - rc(('xtick.major','xtick.minor','ytick.major','ytick.minor'), pad=dist_tick_lab) - - - -def sciy(): - pl.gca().ticklabel_format(style='sci', scilimits=(0,0),axis='y') - -def scix(): - pl.gca().ticklabel_format(style='sci', scilimits=(0,0),axis='x') diff --git a/000_test_round_chamber.py b/tests/000_test_round_chamber.py similarity index 98% rename from 000_test_round_chamber.py rename to tests/000_test_round_chamber.py index 49337cfa..96f84bfb 100644 --- a/000_test_round_chamber.py +++ b/tests/000_test_round_chamber.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW import PyPIC.FiniteDifferences_Staircase_SquareGrid as PIC_FD import PyPIC.FFT_OpenBoundary as PIC_FFT @@ -45,7 +42,7 @@ if PIC_FPPS: picFPPS = PIC_FPPS(200,200,a=R_cham,solverType='Uniform') # build dual grid pic_main = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_main) -pic_dualgrid = AddInternalGrid(pic_main, x_min_internal, x_max_internal, y_min_internal, +pic_dualgrid = AddInternalGrid(pic_main, x_min_internal, x_max_internal, y_min_internal, y_max_internal, Dh_internal, N_nodes_discard) # generate particles diff --git a/002_test_rect_chamber.py b/tests/002_test_rect_chamber.py similarity index 97% rename from 002_test_rect_chamber.py rename to tests/002_test_rect_chamber.py index 256bcbc0..0eea5e97 100644 --- a/002_test_rect_chamber.py +++ b/tests/002_test_rect_chamber.py @@ -1,7 +1,3 @@ -import sys, os -BIN=os.path.expanduser('../') -sys.path.append(BIN) - import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW import PyPIC.FiniteDifferences_Staircase_SquareGrid as PIC_FD import PyPIC.FFT_PEC_Boundary_SquareGrid as PIC_PEC_FFT diff --git a/003_Xmas_tree.py b/tests/003_Xmas_tree.py similarity index 97% rename from 003_Xmas_tree.py rename to tests/003_Xmas_tree.py index 73b79c6b..77f8e2d6 100644 --- a/003_Xmas_tree.py +++ b/tests/003_Xmas_tree.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import pylab as pl import numpy as np from scipy import rand @@ -29,7 +26,7 @@ [1,9], [2,9], [0,11]] - + tree=np.array(tree) x_tree = tree[:,0] y_tree = tree[:,1] @@ -43,9 +40,6 @@ y_tree = np.array([-y_aper]+ list(y_tree)+[y_aper]) - - - x_part = x_aper*(2.*rand(N_part_gen)-1.) y_part = y_aper*(2.*rand(N_part_gen)-1.) @@ -58,14 +52,11 @@ nel_part = 0*x_part+1. - - - chamber = poly.polyg_cham_geom_object({'Vx':na([x_aper, -x_aper, -x_aper, x_aper]), 'Vy':na([y_aper, y_aper, -y_aper, -y_aper]), 'x_sem_ellip_insc':0.99*x_aper, 'y_sem_ellip_insc':0.99*y_aper}) - + picFDSW = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh) picFFTPEC = PIC_PEC_FFT.FFT_PEC_Boundary_SquareGrid(x_aper = chamber.x_aper, y_aper = chamber.y_aper, Dh = Dh) picFFT = PIC_FFT.FFT_OpenBoundary_SquareGrid(x_aper = chamber.x_aper, y_aper = chamber.y_aper, Dh = Dh) @@ -135,7 +126,7 @@ pl.colorbar() pl.savefig('Xmas_efield_open_boudary.png', dpi=200) -if PIC_FPPS: +if PIC_FPPS: x = np.arange(-chamber.x_aper,chamber.x_aper,2*chamber.x_aper/300) y = np.arange(-chamber.y_aper,chamber.y_aper,2*chamber.y_aper/300) X,Y = np.meshgrid(x,y) diff --git a/004_test_gaussian.py b/tests/004_test_gaussian.py similarity index 98% rename from 004_test_gaussian.py rename to tests/004_test_gaussian.py index 292d9392..e9faac4f 100644 --- a/004_test_gaussian.py +++ b/tests/004_test_gaussian.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW import PyPIC.FiniteDifferences_Staircase_SquareGrid as PIC_FD import PyPIC.FFT_OpenBoundary_SquareGrid as PIC_FFT diff --git a/005_testfftw.py b/tests/005_testfftw.py similarity index 95% rename from 005_testfftw.py rename to tests/005_testfftw.py index 3de755e9..07576721 100644 --- a/005_testfftw.py +++ b/tests/005_testfftw.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import pylab as pl import numpy as np from scipy import rand @@ -26,7 +23,7 @@ [1,9], [2,9], [0,11]] - + tree=np.array(tree) x_tree = tree[:,0] y_tree = tree[:,1] @@ -40,9 +37,6 @@ y_tree = np.array([-y_aper]+ list(y_tree)+[y_aper]) - - - x_part = x_aper*(2.*rand(N_part_gen)-1.) y_part = y_aper*(2.*rand(N_part_gen)-1.) @@ -55,15 +49,10 @@ nel_part = 0*x_part+1. - - - chamber = poly.polyg_cham_geom_object({'Vx':na([x_aper, -x_aper, -x_aper, x_aper]), 'Vy':na([y_aper, y_aper, -y_aper, -y_aper]), 'x_sem_ellip_insc':0.99*x_aper, 'y_sem_ellip_insc':0.99*y_aper}) - - picFFT = PIC_FFT.FFT_OpenBoundary_SquareGrid(x_aper = chamber.x_aper, y_aper = chamber.y_aper, Dh = Dh) @@ -79,7 +68,7 @@ for _ in range(N_rep): transf = np.fft.fft2(data) itransf = np.real(np.fft.ifft2(transf*data)) - + t_stop_npfft = time.mktime(time.localtime()) t_npfft = t_stop_npfft-t_start_npfft print('t_npfft', t_npfft) diff --git a/006_time_solve.py b/tests/006_time_solve.py similarity index 95% rename from 006_time_solve.py rename to tests/006_time_solve.py index 7cecad32..3cda774c 100644 --- a/006_time_solve.py +++ b/tests/006_time_solve.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import pylab as pl import numpy as np from scipy import rand @@ -25,7 +22,7 @@ [1,9], [2,9], [0,11]] - + tree=np.array(tree) x_tree = tree[:,0] y_tree = tree[:,1] @@ -39,9 +36,6 @@ y_tree = np.array([-y_aper]+ list(y_tree)+[y_aper]) - - - x_part = x_aper*(2.*rand(N_part_gen)-1.) y_part = y_aper*(2.*rand(N_part_gen)-1.) @@ -54,14 +48,11 @@ nel_part = 0*x_part+1. - - - chamber = poly.polyg_cham_geom_object({'Vx':na([x_aper, -x_aper, -x_aper, x_aper]), 'Vy':na([y_aper, y_aper, -y_aper, -y_aper]), 'x_sem_ellip_insc':0.99*x_aper, 'y_sem_ellip_insc':0.99*y_aper}) - + picFDSW = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh) picFFTPEC = PIC_PEC_FFT.FFT_PEC_Boundary_SquareGrid(x_aper = chamber.x_aper, y_aper = chamber.y_aper, Dh = Dh, fftlib='pyfftw') picFFT = PIC_FFT.FFT_OpenBoundary_SquareGrid(x_aper = chamber.x_aper, y_aper = chamber.y_aper, Dh = Dh, fftlib='pyfftw') diff --git a/007_test_separately.py b/tests/007_test_separately.py similarity index 98% rename from 007_test_separately.py rename to tests/007_test_separately.py index 98650f59..56607ac1 100644 --- a/007_test_separately.py +++ b/tests/007_test_separately.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW import PyPIC.FiniteDifferences_Staircase_SquareGrid as PIC_FD import PyPIC.FFT_OpenBoundary_SquareGrid as PIC_FFT diff --git a/008a_check_interpolation_near_borders.py b/tests/008a_check_interpolation_near_borders.py similarity index 95% rename from 008a_check_interpolation_near_borders.py rename to tests/008a_check_interpolation_near_borders.py index 5cbe83ad..04354af3 100644 --- a/008a_check_interpolation_near_borders.py +++ b/tests/008a_check_interpolation_near_borders.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import numpy as np from PyPIC.geom_impact_ellip import ellip_cham_geom_object @@ -50,12 +47,11 @@ pic.solve(rho = rho_mat) except ValueError as err: print('Got ValueError:', err) - - + Ex_list = [] Ey_list = [] -for xmax_test in xmax_test_list: - +for xmax_test in xmax_test_list: + sc_test = xmax_test/x_aper x_test = sc_test*x_aper*np.cos(thp) y_test = sc_test*y_aper*np.sin(thp) @@ -64,7 +60,7 @@ Ex_list.append(Ex) Ey_list.append(Ey) - + Ex_list = np.array(Ex_list) Ey_list = np.array(Ey_list) diff --git a/008b_check_interpolation_near_borders_plots.py b/tests/008b_check_interpolation_near_borders_plots.py similarity index 97% rename from 008b_check_interpolation_near_borders_plots.py rename to tests/008b_check_interpolation_near_borders_plots.py index d7afe8d9..90c6006d 100644 --- a/008b_check_interpolation_near_borders_plots.py +++ b/tests/008b_check_interpolation_near_borders_plots.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import numpy as np import pylab as pl @@ -28,37 +25,37 @@ erry_rel_list_old = [] erry_rel_list = [] -for ii in range(N_points): +for ii in range(N_points): Ex_ref = ob_ref.Ex[ii,:] Ey_ref = ob_ref.Ey[ii,:] - + Ex_new = ob_new.Ex[ii,:] Ey_new = ob_new.Ey[ii,:] - + Ex_old= ob_old.Ex[ii,:] Ey_old= ob_old.Ey[ii,:] err_abs = np.sqrt(np.sum((Ex_ref-Ex_new)**2+(Ey_ref-Ey_new)**2)) err_rel = err_abs/np.sqrt(np.sum((Ex_ref)**2+(Ey_ref)**2)) - + err_abs_old = np.sqrt(np.sum((Ex_ref-Ex_old)**2+(Ey_ref-Ey_old)**2)) err_rel_old = err_abs_old/np.sqrt(np.sum((Ex_ref)**2+(Ey_ref)**2)) - + erry_rel = np.sqrt(np.sum((Ey_ref-Ey_new)**2))/np.sqrt(np.sum((Ey_ref)**2)) erry_rel_old = np.sqrt(np.sum((Ey_ref-Ey_old)**2))/np.sqrt(np.sum((Ey_ref)**2)) - + err_abs_list.append(err_abs) err_rel_list.append(err_rel) - + err_abs_list_old.append(err_abs_old) err_rel_list_old.append(err_rel_old) - + erry_rel_list_old.append(erry_rel_old) erry_rel_list.append(erry_rel) - + na = np.array -pl.close('all') +pl.close('all') ms.mystyle_arial(fontsz=16, dist_tick_lab=10) pl.figure(1) pl.plot(1000*ob_ref.xmax_test_list, err_abs_list) @@ -181,4 +178,4 @@ pl.savefig(fname+'.png', dpi=200) pl.show() - + diff --git a/009_test_round_chamber_for_new_interp_devel.py b/tests/009_test_round_chamber_for_new_interp_devel.py similarity index 97% rename from 009_test_round_chamber_for_new_interp_devel.py rename to tests/009_test_round_chamber_for_new_interp_devel.py index 4a83521a..b6aa8ae6 100644 --- a/009_test_round_chamber_for_new_interp_devel.py +++ b/tests/009_test_round_chamber_for_new_interp_devel.py @@ -1,6 +1,3 @@ -import sys -sys.path.append('..') - import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW import PyPIC.FFT_OpenBoundary_SquareGrid as PIC_FFT import PyPIC.geom_impact_ellip as ell diff --git a/010_effect_of_tolerances_in_SW.py b/tests/010_effect_of_tolerances_in_SW.py similarity index 92% rename from 010_effect_of_tolerances_in_SW.py rename to tests/010_effect_of_tolerances_in_SW.py index 49ff4aa3..f965d384 100644 --- a/010_effect_of_tolerances_in_SW.py +++ b/tests/010_effect_of_tolerances_in_SW.py @@ -1,9 +1,3 @@ -import sys, os -BIN=os.path.expanduser('../') -sys.path.append(BIN) -BIN=os.path.expanduser('../PyHEADTAIL/testing/script-tests/') -sys.path.append(BIN) - import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW import PyPIC.FiniteDifferences_Staircase_SquareGrid as PIC_FD import PyPIC.geom_impact_ellip as ell @@ -16,14 +10,14 @@ def analytic_solution(x_probes, y_probes, x_part, y_part, nel_part, charge): Ex_probes = [] Ey_probes = [] - + for x_probe, y_probe in zip(x_probes, y_probes): r_probe = np.sqrt(x_probe**2 + y_probe**2) q_inside = np.sum(nel_part[x_part**2+ y_part**2 < r_probe**2])*charge Er_probe = q_inside/(epsilon_0*2*np.pi*r_probe) Ex_probes.append(Er_probe * x_probe/r_probe) Ey_probes.append(Er_probe * y_probe/r_probe) - + return Ex_probes, Ey_probes qe = e @@ -46,7 +40,7 @@ def analytic_solution(x_probes, y_probes, x_part, y_part, nel_part, charge): chamber = ell.ellip_cham_geom_object(x_aper = x_aper, y_aper = y_aper) # build main pic -pic_SW = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, +pic_SW = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_main, tol_stem = tol_stem, tol_der = tol_der) @@ -100,5 +94,3 @@ def analytic_solution(x_probes, y_probes, x_part, y_part, nel_part, charge): pl.show() - - diff --git a/011_LHC_multigrid_test.py b/tests/011_LHC_multigrid_test.py similarity index 90% rename from 011_LHC_multigrid_test.py rename to tests/011_LHC_multigrid_test.py index 75c84f6f..ee5a6043 100644 --- a/011_LHC_multigrid_test.py +++ b/tests/011_LHC_multigrid_test.py @@ -1,8 +1,3 @@ -import sys, os -BIN=os.path.expanduser('../') -sys.path.append(BIN) -BIN=os.path.expanduser('../PyHEADTAIL/testing/script-tests/') -sys.path.append(BIN) from LHC import LHC import PyPIC.geom_impact_ellip as ell import PyPIC.FiniteDifferences_ShortleyWeller_SquareGrid as PIC_FDSW @@ -34,7 +29,7 @@ # generate beam -bunch = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = intensity, +bunch = machine.generate_6D_Gaussian_bunch(n_macroparticles = n_macroparticles, intensity = intensity, epsn_x = epsn_x, epsn_y = epsn_y, sigma_z = sigma_z) @@ -59,9 +54,9 @@ Vx, Vy = chamber.points_on_boundary(N_points=200) # build Bassetti Erskine -pic_BE = PIC_BE.Interpolated_Bassetti_Erskine(x_aper=x_aper, y_aper=y_aper, Dh=Dh_BE, sigmax=bunch.sigma_x(), sigmay=bunch.sigma_y(), +pic_BE = PIC_BE.Interpolated_Bassetti_Erskine(x_aper=x_aper, y_aper=y_aper, Dh=Dh_BE, sigmax=bunch.sigma_x(), sigmay=bunch.sigma_y(), n_imag_ellip=20, tot_charge=bunch.intensity*bunch.charge) - + # build single grid pic pic_singlegrid = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_single, sparse_solver = sparse_solver) @@ -69,13 +64,13 @@ pic_singlegrid_ext = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_single_ext, sparse_solver = sparse_solver) # build telescope -pic_multigrid = AddTelescopicGrids(pic_main = pic_singlegrid_ext, f_telescope = 0.3, - target_grid = {'x_min_target':-Sx_target/2., 'x_max_target':Sx_target/2.,'y_min_target':-Sy_target/2.,'y_max_target':Sy_target/2.,'Dh_target':Dh_target}, +pic_multigrid = AddTelescopicGrids(pic_main = pic_singlegrid_ext, f_telescope = 0.3, + target_grid = {'x_min_target':-Sx_target/2., 'x_max_target':Sx_target/2.,'y_min_target':-Sy_target/2.,'y_max_target':Sy_target/2.,'Dh_target':Dh_target}, N_nodes_discard = 3., N_min_Dh_main = 10, sparse_solver=sparse_solver) pic_singlegrid.scatter(bunch.x, bunch.y, bunch.particlenumber_per_mp+bunch.y*0., charge=qe) -pic_multigrid.scatter(bunch.x, bunch.y, bunch.particlenumber_per_mp+bunch.y*0., charge=qe) - +pic_multigrid.scatter(bunch.x, bunch.y, bunch.particlenumber_per_mp+bunch.y*0., charge=qe) + #scatter and solve #pic solve timing import time @@ -85,10 +80,10 @@ for _ in range(N_rep_test_single): pic_singlegrid.solve() t_stop = time.mktime(time.localtime()) -t_sw_single = (t_stop-t_start)/N_rep_test_single +t_sw_single = (t_stop-t_start)/N_rep_test_single print('solving time singlegrid ', t_sw_single) -N_rep_test_multi = 10000 +N_rep_test_multi = 10000 print('Solving PIC multi %d times'%N_rep_test_multi) t_start = time.mktime(time.localtime()) for _ in range(N_rep_test_multi): @@ -102,7 +97,7 @@ theta=np.linspace(0., 2*np.pi, 1000) n_sigma_probes = 1. x_probes = n_sigma_probes*bunch.sigma_x()*np.cos(theta) -y_probes = n_sigma_probes*bunch.sigma_y()*np.sin(theta) +y_probes = n_sigma_probes*bunch.sigma_y()*np.sin(theta) # get field at probes Ex_BE, Ey_BE = pic_BE.gather(x_probes, y_probes) @@ -126,7 +121,7 @@ pl.plot(x_probes, y_probes, 'c--', label = 'probe') pl.xlabel('x [m]') pl.ylabel('y [m]') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') +pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') pl.axis('equal') pl.legend(loc='best') @@ -136,7 +131,7 @@ pl.plot(theta*180/np.pi, Ex_multigrid, '.-r', label = 'Multigrid') pl.xlabel('theta[deg]') pl.ylabel('Ex [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') +pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') pl.grid() pl.legend(loc='best') @@ -146,7 +141,7 @@ pl.plot(theta*180/np.pi, Ey_multigrid, '.-r', label = 'Multigrid') pl.xlabel('theta[deg]') pl.ylabel('Ey [V/m] ') -pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') +pl.ticklabel_format(style='sci', scilimits=(0,0),axis='x') pl.ticklabel_format(style='sci', scilimits=(0,0),axis='y') pl.grid() pl.legend(loc='best') @@ -177,7 +172,7 @@ pl.loglog(r_probes_val, RMSE_multigrid, '.-g', label = 'Multigrid vs BE ') pl.xlabel('r [m]') pl.ylabel('RMS error') -pl.title('$\sigma_x$ = %.2e [m]\n $\sigma_y$ = %.2e [m] \n $\Delta h_{single}$ = %.2e [m]\n $\Delta h_{multi}$ = %.2e [m]\n $\Delta h_{BE}$ = %.2e [m]\n Solving time: $t_{single}$ = %.1f ms, $t_{multi}$ = %.1f ms'%(bunch.sigma_x(), bunch.sigma_y(), Dh_single, +pl.title('$\sigma_x$ = %.2e [m]\n $\sigma_y$ = %.2e [m] \n $\Delta h_{single}$ = %.2e [m]\n $\Delta h_{multi}$ = %.2e [m]\n $\Delta h_{BE}$ = %.2e [m]\n Solving time: $t_{single}$ = %.1f ms, $t_{multi}$ = %.1f ms'%(bunch.sigma_x(), bunch.sigma_y(), Dh_single, Dh_target, Dh_BE, t_sw_single*1000., t_sw_multi*1000.)) pl.subplots_adjust(bottom = .13, top = .70) pl.grid() @@ -207,7 +202,7 @@ pl.loglog(n_sigma_probes, RMSE_multigrid, '.-g', label = 'Multigrid vs BE ') pl.xlabel('$\sigma$') pl.ylabel('RMS error') -pl.title('$\sigma_x$ = %.2e [m]\n $\sigma_y$ = %.2e [m] \n $\Delta h_{single}$ = %.2e [m]\n $\Delta h_{multi}$ = %.2e [m]\n $\Delta h_{BE}$ = %.2e [m]\n Solving time: $t_{single}$ = %.1f ms, $t_{multi}$ = %.1f ms'%(bunch.sigma_x(), bunch.sigma_y(), Dh_single, +pl.title('$\sigma_x$ = %.2e [m]\n $\sigma_y$ = %.2e [m] \n $\Delta h_{single}$ = %.2e [m]\n $\Delta h_{multi}$ = %.2e [m]\n $\Delta h_{BE}$ = %.2e [m]\n Solving time: $t_{single}$ = %.1f ms, $t_{multi}$ = %.1f ms'%(bunch.sigma_x(), bunch.sigma_y(), Dh_single, Dh_target, Dh_BE, t_sw_single*1000., t_sw_multi*1000.)) pl.subplots_adjust(bottom = .13, top = .70) pl.grid() @@ -224,21 +219,21 @@ yn=yn.T yn=yn.flatten() #pic gather -Ex_BE_n, Ey_BE_n = pic_BE.gather(xn, yn) +Ex_BE_n, Ey_BE_n = pic_BE.gather(xn, yn) Ex_BE_matrix=np.reshape(Ex_BE_n,(len(y_grid_probes),len(x_grid_probes)), 'F').T Ey_BE_matrix=np.reshape(Ey_BE_n,(len(y_grid_probes),len(x_grid_probes)), 'F').T -Ex_singlegrid_n, Ey_singlegrid_n = pic_singlegrid.gather(xn, yn) +Ex_singlegrid_n, Ey_singlegrid_n = pic_singlegrid.gather(xn, yn) Ex_singlegrid_matrix=np.reshape(Ex_singlegrid_n,(len(y_grid_probes),len(x_grid_probes)), 'F').T Ey_singlegrid_matrix=np.reshape(Ey_singlegrid_n,(len(y_grid_probes),len(x_grid_probes)), 'F').T -Ex_multigrid_n, Ey_multigrid_n = pic_multigrid.gather(xn, yn) +Ex_multigrid_n, Ey_multigrid_n = pic_multigrid.gather(xn, yn) Ex_multigrid_matrix=np.reshape(Ex_multigrid_n,(len(y_grid_probes),len(x_grid_probes)), 'F').T Ey_multigrid_matrix=np.reshape(Ey_multigrid_n,(len(y_grid_probes),len(x_grid_probes)), 'F').T pl.figure(4, figsize=(12, 6)).patch.set_facecolor('w') sp1 = pl.subplot(121) -pl.pcolormesh(x_grid_probes, y_grid_probes, +pl.pcolormesh(x_grid_probes, y_grid_probes, np.log10(np.sqrt((((Ex_singlegrid_matrix-Ex_BE_matrix)**2+(Ey_singlegrid_matrix-Ey_BE_matrix)**2)/(Ex_BE_matrix**2+Ey_BE_matrix**2)))).T, vmax=0., vmin=-7.0) pl.title('RMS error Singlegrid - BE') @@ -248,11 +243,11 @@ cb.formatter.set_powerlimits((0, 0)) cb.update_ticks() cb.set_label('RMS error') -sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='x') +sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='x') sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='y') sp2 = pl.subplot(122, sharex = sp1, sharey = sp1) -pl.pcolormesh(x_grid_probes, y_grid_probes, +pl.pcolormesh(x_grid_probes, y_grid_probes, np.log10(np.sqrt((((Ex_multigrid_matrix-Ex_BE_matrix)**2+(Ey_multigrid_matrix-Ey_BE_matrix)**2)/(Ex_BE_matrix**2+Ey_BE_matrix**2)))).T, vmax=0., vmin=-7.0) pl.title('RMS error Multigrid - BE') @@ -262,7 +257,7 @@ cb.formatter.set_powerlimits((0, 0)) cb.update_ticks() cb.set_label('RMS error') -sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='x') +sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='x') sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='y') pl.subplots_adjust(bottom = .13,top = .70) pl.suptitle('$\sigma_x$ = %.2e [m]\n $\sigma_y$ = %.2e [m] \n $\Delta h_{single}$ = %.2e [m]\n $\Delta h_{multi}$ = %.2e [m]\n $\Delta h_{BE}$ = %.2e [m]\n Solving time: $t_{single}$ = %.1f ms, $t_{multi}$ = %.1f ms'%(bunch.sigma_x(), bunch.sigma_y(), Dh_single, diff --git a/012_test_states.py b/tests/012_test_states.py similarity index 94% rename from 012_test_states.py rename to tests/012_test_states.py index ea3385c9..883b118d 100644 --- a/012_test_states.py +++ b/tests/012_test_states.py @@ -1,7 +1,3 @@ -import sys -sys.path.append('../') - - import numpy as np import pylab as pl from scipy.constants import e as qe @@ -46,8 +42,8 @@ Dh_single = .5e-3 pic_singlegrid = PIC_FDSW.FiniteDifferences_ShortleyWeller_SquareGrid(chamb = chamber, Dh = Dh_single, sparse_solver = 'PyKLU') from PyPIC.MultiGrid import AddTelescopicGrids -pic_list.append(AddTelescopicGrids(pic_main = pic_singlegrid, f_telescope = 0.3, - target_grid = {'x_min_target':-Sx_target/2., 'x_max_target':Sx_target/2.,'y_min_target':-Sy_target/2.,'y_max_target':Sy_target/2.,'Dh_target':Dh_target}, +pic_list.append(AddTelescopicGrids(pic_main = pic_singlegrid, f_telescope = 0.3, + target_grid = {'x_min_target':-Sx_target/2., 'x_max_target':Sx_target/2.,'y_min_target':-Sy_target/2.,'y_max_target':Sy_target/2.,'Dh_target':Dh_target}, N_nodes_discard = 3., N_min_Dh_main = 10, sparse_solver='PyKLU')) # test: @@ -56,7 +52,7 @@ ms.mystyle_arial(fontsz = 14) for i_pic, pic in enumerate(pic_list): - + # standard solve and gather pic.scatter(x_mp, y_mp, nel_mp, charge = qe) pic.solve() @@ -69,10 +65,10 @@ list_states = [] for _ in fact_states: list_states.append(pic.get_state_object()) - + for i_state, state in enumerate(list_states): state.scatter(x_mp, y_mp, nel_mp*fact_states[i_state], charge = qe) - + # solve states pic.solve_states(list_states) @@ -109,11 +105,11 @@ pl.subplot(2,1,2) pl.plot(theta_probes, Ey_prb_single_state, '.', color=colorcurr, label = 'Single state') pl.plot(theta_probes, Ey_probes*fact_single_state, '-', color=colorcurr, label = 'Single ref.') - - sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='x') + + sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='x') sp1.ticklabel_format(style='sci', scilimits=(0,0),axis='y') sp1.legend(loc='center left', bbox_to_anchor=(1, 0.)) - sp2.ticklabel_format(style='sci', scilimits=(0,0),axis='x') + sp2.ticklabel_format(style='sci', scilimits=(0,0),axis='x') sp2.ticklabel_format(style='sci', scilimits=(0,0),axis='y') pl.subplots_adjust(left = .10, right = .77, bottom = .13,top = .90, hspace = .40) pl.suptitle(str(pic.__class__).replace('<','').replace('>','').split('.')[-1].replace("'", '')) diff --git a/tests/LHC.py b/tests/LHC.py new file mode 100644 index 00000000..063223f0 --- /dev/null +++ b/tests/LHC.py @@ -0,0 +1,190 @@ +import numpy as np +from scipy.constants import c, e, m_p + +from PyHEADTAIL.machines.synchrotron import Synchrotron + + +class EmptyObject(object): + pass + + +class LHC(Synchrotron): + + def __init__(self, machine_configuration=None, optics_mode='smooth', **kwargs): + + + pp = EmptyObject() + pp.machine_configuration = machine_configuration + pp.optics_mode = optics_mode + + pp.longitudinal_mode = 'non-linear' + pp.alpha = 3.225e-04 + pp.h_RF = 35640 + pp.mass = m_p + pp.charge = e + pp.RF_at = 'middle' + + if pp.machine_configuration == 'Injection': + pp.p0 = 450e9 * e / c + pp.p_increment = 0. + pp.accQ_x = 64.28 + pp.accQ_y = 59.31 + pp.V_RF = 6e6 + pp.dphi_RF = 0. + elif machine_configuration == '6.5_TeV_collision_tunes': + pp.p0 = 6500e9 * e / c + pp.p_increment = 0. + pp.accQ_x = 64.31 + pp.accQ_y = 59.32 + pp.V_RF = 12e6 + pp.dphi_RF = 0. + else: + raise ValueError('machine_configuration not recognized!') + + if pp.optics_mode == 'smooth': + if 's' in list(kwargs.keys()): + raise ValueError('s vector cannot be provided if optics_mode = "smooth"') + + pp.n_segments = kwargs['n_segments'] + pp.circumference = 26658.8832 + + pp.name = None + + pp.beta_x = 92.7 + pp.D_x = 0 + pp.beta_y = 93.2 + pp.D_y = 0 + + pp.alpha_x = None + pp.alpha_y = None + + pp.s = None + + elif pp.optics_mode == 'non-smooth': + if 'n_segments' in list(kwargs.keys()): + raise ValueError('n_segments cannot be provided if optics_mode = "non-smooth"') + pp.n_segments = None + pp.circumference = None + + pp.name = kwargs['name'] + + pp.beta_x = kwargs['beta_x'] + pp.beta_y = kwargs['beta_y'] + + try: + pp.D_x = kwargs['D_x'] + except KeyError: + pp.D_x = 0 * np.array(kwargs['s']) + try: + pp.D_y = kwargs['D_y'] + except KeyError: + pp.D_y = 0 * np.array(kwargs['s']) + + pp.alpha_x = kwargs['alpha_x'] + pp.alpha_y = kwargs['alpha_y'] + + pp.s = kwargs['s'] + + else: + raise ValueError('optics_mode not recognized!') + + # detunings + pp.Qp_x = 0 + pp.Qp_y = 0 + + pp.app_x = 0 + pp.app_y = 0 + pp.app_xy = 0 + + pp.i_octupole_focusing = None + pp.i_octupole_defocusing = None + pp.octupole_knob = None + + for attr in list(kwargs.keys()): + if kwargs[attr] is not None: + if type(kwargs[attr]) is list or type(kwargs[attr]) is np.ndarray: + str2print = '[%s ...]'%repr(kwargs[attr][0]) + else: + str2print = repr(kwargs[attr]) + self.prints('Synchrotron init. From kwargs: %s = %s' + % (attr, str2print)) + if not hasattr(pp, attr): + raise NameError("I don't understand %s"%attr) + + setattr(pp, attr, kwargs[attr]) + + + + if pp.i_octupole_focusing is not None or pp.i_octupole_defocusing is not None: + if pp.octupole_knob is not None: + raise ValueError('octupole_knobs and octupole currents cannot be used at the same time!') + pp.app_x, pp.app_y, pp.app_xy = self._anharmonicities_from_octupole_current_settings( + pp.i_octupole_focusing, pp.i_octupole_defocusing) + self.i_octupole_focusing = pp.i_octupole_focusing + self.i_octupole_defocusing = pp.i_octupole_defocusing + + if pp.octupole_knob is not None: + if pp.i_octupole_focusing is not None or pp.i_octupole_defocusing is not None: + raise ValueError('octupole_knobs and octupole currents cannot be used at the same time!') + pp.i_octupole_focusing, pp.i_octupole_defocusing = self._octupole_currents_from_octupole_knobs(pp.octupole_knob, pp.p0) + pp.app_x, pp.app_y, pp.app_xy = self._anharmonicities_from_octupole_current_settings( + pp.i_octupole_focusing, pp.i_octupole_defocusing) + self.i_octupole_focusing = pp.i_octupole_focusing + self.i_octupole_defocusing = pp.i_octupole_defocusing + + + super(LHC, self).__init__(optics_mode=pp.optics_mode, circumference=pp.circumference, n_segments=pp.n_segments, + s=pp.s, name=pp.name, + alpha_x=pp.alpha_x, beta_x=pp.beta_x, D_x=pp.D_x, alpha_y=pp.alpha_y, beta_y=pp.beta_y, D_y=pp.D_y, + accQ_x=pp.accQ_x, accQ_y=pp.accQ_y, Qp_x=pp.Qp_x, Qp_y=pp.Qp_y, app_x=pp.app_x, app_y=pp.app_y, app_xy=pp.app_xy, + alpha_mom_compaction=pp.alpha, longitudinal_mode=pp.longitudinal_mode, + h_RF=np.atleast_1d(pp.h_RF), V_RF=np.atleast_1d(pp.V_RF), dphi_RF=np.atleast_1d(pp.dphi_RF), + p0=pp.p0, p_increment=pp.p_increment, + charge=pp.charge, mass=pp.mass, RF_at=pp.RF_at) + + + def _anharmonicities_from_octupole_current_settings(self, i_octupole_focusing, i_octupole_defocusing): + """Calculate the constants of proportionality app_x, app_y and + app_xy (== app_yx) for the amplitude detuning introduced by the + LHC octupole magnets (aka. LHC Landau octupoles) from the + electric currents i_octupole_focusing [A] and i_octupole_defocusing [A] flowing + through the magnets. The maximum current is given by + i_max = +/- 550 [A]. The values app_x, app_y, app_xy obtained + from the formulae are proportional to the strength of detuning + for one complete turn around the accelerator, i.e. one-turn + values. + The calculation is based on formulae (3.6) taken from 'The LHC + transverse coupled-bunch instability' by N. Mounet, EPFL PhD + Thesis, 2012. Values (hard-coded numbers below) are valid for + LHC Landau octupoles before LS1. Beta functions in x and y are + correctly taken into account. Note that here, the values of + app_x, app_y and app_xy are not normalized to the reference + momentum p0. This is done only during the calculation of the + detuning in the corresponding detune method of the + AmplitudeDetuningSegment. + More detailed explanations and references on how the formulae + were obtained are given in the PhD thesis (pg. 85ff) cited + above. + """ + i_max = 550. # [A] + E_max = 7000. # [GeV] + + app_x = E_max * (267065. * i_octupole_focusing / i_max - + 7856. * i_octupole_defocusing / i_max) + app_y = E_max * (9789. * i_octupole_focusing / i_max - + 277203. * i_octupole_defocusing / i_max) + app_xy = E_max * (-102261. * i_octupole_focusing / i_max + + 93331. * i_octupole_defocusing / i_max) + + # Convert to SI units. + convert_to_SI = e / (1.e-9 * c) + app_x *= convert_to_SI + app_y *= convert_to_SI + app_xy *= convert_to_SI + + return app_x, app_y, app_xy + + def _octupole_currents_from_octupole_knobs(self, octupole_knob, p0): + i_octupole_focusing = 19.557 * octupole_knob / (-1.5) * p0 / 2.4049285931335872e-16 + i_octupole_defocusing = - i_octupole_focusing + return i_octupole_focusing, i_octupole_defocusing