From 921ca9fc1f89deb05cbf5e70c29384a7906a42a2 Mon Sep 17 00:00:00 2001 From: David Sagan Date: Tue, 24 Oct 2023 21:18:00 -0400 Subject: [PATCH] Better set_tune_via_group_knobs algorithm. --- bmad/code/set_tune_via_group_knobs.f90 | 23 ++++-- bmad/multiparticle/beam_file_io.f90 | 99 ++++++++++++++++++++++++-- bsim/modules/ts_mod.f90 | 38 +++++++--- tao/version/tao_version_mod.f90 | 2 +- 4 files changed, 144 insertions(+), 18 deletions(-) diff --git a/bmad/code/set_tune_via_group_knobs.f90 b/bmad/code/set_tune_via_group_knobs.f90 index 5c3db7eb89..7be50dc0a3 100644 --- a/bmad/code/set_tune_via_group_knobs.f90 +++ b/bmad/code/set_tune_via_group_knobs.f90 @@ -31,7 +31,7 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err) real(rp) phi_set(2), dphi_a, dphi_b, dQ_max real(rp) phi_a, phi_b, d_a1, d_a2, d_b1, d_b2, det -real(rp) d1, d2, del0 +real(rp) d1, d2, del0, factor real(rp) :: phi_array(2) real(rp), allocatable :: deriv1(:), deriv2(:), kinit(:) @@ -48,6 +48,7 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err) dQ_max = 0.001 del0 = 0.001 ok = .false. +factor = 1 rf_on = rf_is_on(branch) allocate(kinit(branch%n_ele_track), deriv1(branch%n_ele_track), deriv2(branch%n_ele_track)) @@ -58,7 +59,7 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err) ! Q tune -do i = 1, 10 +do i = 1, 20 call lattice_bookkeeper(branch%lat) if (rf_on) then @@ -72,7 +73,19 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err) call lat_make_mat6 (branch%lat, -1, orb, branch%ix_branch) call twiss_at_start(branch%lat, status, branch%ix_branch, print_err) - if (status /= ok$) return + if (status /= ok$) then + if (i == 1) return + group1%control%var(1)%value = group1%control%var(1)%value - factor * d1 + group2%control%var(1)%value = group2%control%var(1)%value - factor * d2 + factor = 0.5 * factor + group1%control%var(1)%value = group1%control%var(1)%value + factor * d1 + group2%control%var(1)%value = group2%control%var(1)%value + factor * d2 + call set_flags_for_changed_attribute (group1) + call set_flags_for_changed_attribute (group2) + cycle + else + factor = 1.0 + endif call twiss_propagate_all (branch%lat, branch%ix_branch, err) if (err) return @@ -113,8 +126,8 @@ function set_tune_via_group_knobs (phi_set, branch, group_knobs, orb, print_err) ! Put in the changes - group1%control%var(1)%value = group1%control%var(1)%value + d1 - group2%control%var(1)%value = group2%control%var(1)%value + d2 + group1%control%var(1)%value = group1%control%var(1)%value + factor * d1 + group2%control%var(1)%value = group2%control%var(1)%value + factor * d2 call set_flags_for_changed_attribute (group1) call set_flags_for_changed_attribute (group2) diff --git a/bmad/multiparticle/beam_file_io.f90 b/bmad/multiparticle/beam_file_io.f90 index 139945c39f..6639210105 100644 --- a/bmad/multiparticle/beam_file_io.f90 +++ b/bmad/multiparticle/beam_file_io.f90 @@ -10,22 +10,26 @@ module beam_file_io !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !+ -! Subroutine write_beam_file (file_name, beam, new_file, file_format, lat) +! Subroutine write_beam_file (file_name, beam, new_file, file_format, lat, cols) ! ! Routine to write a beam file. ! ! A '.h5' suffix will be appended to the created file if hdf5$ format is used and file_name does not ! already have a '.h5' or '.hdf5' suffix. ! +! The ASCII format is ASCII::4 except if cols is present. In this case ASCII::5 is used. +! ! Input: ! file_name -- character(*): Name of file. ! beam -- beam_struct: Beam to write ! new_file -- logical, optional: New file or append? Default = True. ! file_format -- logical, optional: ascii$, or hdf5$ (default). ! lat -- lat_struct, optional: If present, lattice info will be writen to hdf5 files. +! cols -- character(*): Table columns. Possibilities are coord_struct component names. +! column names should be separated by spaces. !- -subroutine write_beam_file (file_name, beam, new_file, file_format, lat) +subroutine write_beam_file (file_name, beam, new_file, file_format, lat, cols) type (beam_struct), target :: beam type (bunch_struct), pointer :: bunch @@ -36,6 +40,7 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat) integer, optional :: file_format character(*) file_name +character(*), optional :: cols character(200) full_name character(*), parameter :: r_name = 'write_beam_file' @@ -58,6 +63,11 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat) ! +if (present(cols)) then + call write_ascii4_beam_file(file_name, beam, new_file, cols, lat) + return +endif + if (logic_option(.true., new_file)) then open (iu, file = full_name) write (iu, '(a)') '!ASCII::3' @@ -87,6 +97,64 @@ subroutine write_beam_file (file_name, beam, new_file, file_format, lat) end subroutine write_beam_file +!----------------------------------------------------------------------------- +!----------------------------------------------------------------------------- +!----------------------------------------------------------------------------- +!+ +! Subroutine write_ascii4_beam_file (file_name, beam, new_file, cols, lat) +! +! Routine to write a beam file in ASCII::4 format. +! +! A '.h5' suffix will be appended to the created file if hdf5$ format is used and file_name does not +! already have a '.h5' or '.hdf5' suffix. +! +! The ASCII format is ASCII::4 except if cols is present. In this case ASCII::5 is used. +! +! Input: +! file_name -- character(*): Name of file. +! beam -- beam_struct: Beam to write +! new_file -- logical, optional: New file or append? Default = True. +! cols -- character(*): Table columns. Possibilities are coord_struct component names. +! column names should be separated by spaces. +! lat -- lat_struct, optional: If present, lattice info will be writen to hdf5 files. +!- + +subroutine write_ascii4_beam_file (file_name, beam, new_file, cols, lat) + +type (beam_struct), target :: beam +type (bunch_struct), pointer :: bunch +type (coord_struct), pointer :: p +type (lat_struct), optional :: lat + +integer j, iu, ib, ip, ix_ele, n, n0 + +character(*) file_name +character(*), optional :: cols +character(200) full_name +character(*), parameter :: r_name = 'write_ascii4_beam_file' + +logical, optional :: new_file +logical error, append + +! + +iu = lunget() +call fullfilename (file_name, full_name) + +if (logic_option(.true., new_file)) then + open (iu, file = full_name) + write (iu, '(a)') '!ASCII::3' +else + open (iu, file = full_name, access = 'append') +endif + +! + + +close (iu) + +end subroutine write_ascii4_beam_file + !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- @@ -636,10 +704,15 @@ subroutine read_beam_ascii4 (iu, file_name, beam, beam_init, err_flag, ele, prin call read_switch(line(:ix), p0%species, str, err) if (err) return - case ('s_position'); p0%s = read_param(line) - case ('time'); p0%t = read_param(line) + case ('r'); p0%r = read_param(line) + case ('s'); p0%s = read_param(line) + case ('t'); p0%t = read_param(line) case ('p0c'); p0%p0c = read_param(line) + case ('dt_ref'); p0%dt_ref = read_param(line) case ('charge'); p0%charge = read_param(line) + case ('spin'); call read_params(line, p0%spin) + case ('field'); call read_params(line, p0%field) + case ('phase'); call read_params(line, p0%phase) case ('time_dir'); p0%time_dir = nint(read_param(line)) case ('direction'); p0%direction = nint(read_param(line)) case ('ix_ele'); p0%ix_ele = nint(read_param(line)) @@ -839,6 +912,24 @@ end function read_param !--------------------------------------------------------------------------------------------------- ! contains +subroutine read_params(line, param) +character(*) line +real(rp) param(:) +integer ix, ios + +! + +ix = index(line, '=') +read(line(ix+1:), *, iostat = ios) param +if (ios /= 0 .or. ix == 0) then + call out_io (s_error$, r_name, 'ERROR READING BEAM FILE PARAMETER!') +endif + +end subroutine read_params + +!--------------------------------------------------------------------------------------------------- +! contains + function read_string(line) result (str) character(*) line character(200) str diff --git a/bsim/modules/ts_mod.f90 b/bsim/modules/ts_mod.f90 index ffeb475990..30b0df3fe7 100644 --- a/bsim/modules/ts_mod.f90 +++ b/bsim/modules/ts_mod.f90 @@ -14,11 +14,11 @@ module ts_mod character(100) :: lat_file = '', dat_out_file = '', quad_mask = '' character(40) :: group_knobs(2) = ["", ""] character(40) :: test = "" - real(rp) :: Q_a0 = 0, Q_a1 = 0, dQ_a = 0 - real(rp) :: Q_b0 = 0, Q_b1 = 0, dQ_b = 0 - real(rp) :: Q_z0 = 0, Q_z1 = 0, dQ_z = 0 + real(rp) :: Q_a0 = real_garbage$, Q_a1 = real_garbage$, dQ_a = real_garbage$ + real(rp) :: Q_b0 = real_garbage$, Q_b1 = real_garbage$, dQ_b = real_garbage$ + real(rp) :: Q_z0 = real_garbage$, Q_z1 = real_garbage$, dQ_z = real_garbage$ real(rp) :: pz0 = 0, pz1 = 0, dpz = 0 - real(rp) :: a_emit = 0, b_emit = 0, sig_pz + real(rp) :: a_emit = 0, b_emit = 0, sig_pz = 0 real(rp) :: a0_amp = 0, b0_amp = 0, pz0_amp = 0 real(rp) :: timer_print_dtime = 120 integer :: n_turn = 0 @@ -80,13 +80,31 @@ subroutine ts_init_params (ts, ts_com) ts_com%master_input_file = 'tune_scan.init' if (n_arg == 1) call get_command_argument(1, ts_com%master_input_file) +bmad_com%auto_bookkeeper = .false. + open (unit= 1, file = ts_com%master_input_file, status = 'old', action = 'read') read(1, nml = params) close (1) -ts%Q_z0 = abs(ts%Q_z0) -ts%Q_z1 = abs(ts%Q_z1) -ts%dQ_z = abs(ts%dQ_z) +if (ts%Q_a0 == real_garbage$ .or. ts%Q_a1 == real_garbage$ .or. ts%dQ_a == real_garbage$) then + print '(a)', 'Error: One of ts%Q_a0, ts%Q_a1, ts%dQ_a is not set! Stopping here.' + stop +endif + +if (ts%Q_b0 == real_garbage$ .or. ts%Q_b1 == real_garbage$ .or. ts%dQ_b == real_garbage$) then + print '(a)', 'Error: One of ts%Q_b0, ts%Q_b1, ts%dQ_b is not set! Stopping here.' + stop +endif + +if (ts%rf_on) then + if (ts%Q_z0 == real_garbage$ .or. ts%Q_z1 == real_garbage$ .or. ts%dQ_z == real_garbage$) then + print '(a)', 'Error: One of ts%Q_z0, ts%Q_z1, ts%dQ_z is not set! Stopping here.' + stop + endif + ts%Q_z0 = abs(ts%Q_z0) + ts%Q_z1 = abs(ts%Q_z1) + ts%dQ_z = abs(ts%dQ_z) +endif if (ts%dat_out_file == '') call file_suffixer(ts_com%master_input_file, ts%dat_out_file, 'dat', .true.) @@ -201,7 +219,7 @@ recursive subroutine ts_track_particle(ts, ts_com, jtune, ts_dat) type (ele_struct), pointer :: ele real(rp) Jvec0(1:6), Jvec(1:6), init_vec(6), amp(3), r(3), v_mat(4,4) -integer jtune(3), nt, track_state, status +integer jtune(3), nt, track_state, status, iqm logical ok, error ! @@ -213,8 +231,10 @@ recursive subroutine ts_track_particle(ts, ts_com, jtune, ts_dat) ts_dat%tune(2) = ts_com%int_Qb + ts%Q_b0 + jtune(2)*ts%dQ_b if (ts%rf_on) then ts_dat%tune(3) = -(ts%Q_z0 + jtune(3)*ts%dQ_z) + iqm = 3 else ts_dat%tune(3) = ts%pz0 + jtune(3)*ts%dpz + iqm = 2 endif ring = ts_com%ring ! Use copy in case tune setting fails, which may garble the lattice @@ -222,6 +242,8 @@ recursive subroutine ts_track_particle(ts, ts_com, jtune, ts_dat) ele => ring%ele(0) ok = set_tune_3d (ring%branch(0), ts_dat%tune, ts%quad_mask, ts%use_phase_trombone, ts%rf_on, ts%group_knobs) ! Tunes in radians. +if (.not. ok) print '(a, 3f9.4)', 'Cannot set tunes (due to resonance?) for: ', ts_dat%tune(1:iqm) + if (.not. ok) return ! Tunes could not be set, probably on a resonance. if (ts%rf_on) then diff --git a/tao/version/tao_version_mod.f90 b/tao/version/tao_version_mod.f90 index 97e70d83a2..2ea94a5a3a 100644 --- a/tao/version/tao_version_mod.f90 +++ b/tao/version/tao_version_mod.f90 @@ -6,5 +6,5 @@ !- module tao_version_mod -character(*), parameter :: tao_version_date = "2023/10/16 23:05:40" +character(*), parameter :: tao_version_date = "2023/10/24 11:25:06" end module