From 23398eca869bf2b9b104cdc8d5aec3360a3940b6 Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 20 Sep 2023 15:25:41 +0200 Subject: [PATCH 1/8] add constant detuning term --- pywit/landau_damping.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/pywit/landau_damping.py b/pywit/landau_damping.py index 97127573..f91dc0ab 100644 --- a/pywit/landau_damping.py +++ b/pywit/landau_damping.py @@ -60,6 +60,7 @@ def dispersion_integral_2d(tune_shift: np.ndarray, b_direct: float, b_cross: flo def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, b_direct_ref: float, b_cross_ref: float, + b_direct_add: float = 0, b_cross_add: float = 0, fraction_of_qs_allowed_on_positive_side: float = 0.05, distribution: str = 'gaussian', tolerance=1e-10): """ @@ -72,6 +73,10 @@ def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, b_direct_ref the x plane or $\alpha_y \sigma_y$ if working in the y plane) :param b_cross_ref: the cross detuning coefficient multiplied by sigma (i.e. $\alpha_{xy} \sigma_y$ if working in the x plane or $\alpha_{yx} \sigma_x$ if working in the y plane) + :param b_direct_add: an additional contribution to the direct detuning coefficient which is kept fixed in the + procedure to find the stability threshold + :param b_cross_add: an additional contribution to the cross detuning coefficient which is kept fixed in the + procedure to find the stability threshold :param distribution: the transverse distribution of the beam. It can be 'gaussian' or 'parabolic' :param fraction_of_qs_allowed_on_positive_side: to determine azimuthal mode number l_mode (around which is drawn the stability diagram), one can consider positive tune shift up to this fraction of q_s (default=5%) @@ -93,8 +98,8 @@ def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, b_direct_ref # function to solve (distance in imag. part w.r.t stab. diagram, as a function of oct. current) def f(b_direct): - b_direct_i = b_direct - b_cross_i = b_ratio * b_direct + b_direct_i = b_direct + b_direct_add + b_cross_i = b_ratio * b_direct + b_cross_add stab = [-1. / dispersion_integral_2d(t_s, b_direct_i, b_cross_i, distribution=distribution) for e in (-1, 1) for t_s in b_direct_i * e * 10. ** np.arange(-3, 2, 0.01)[::e]] # note: one has to reverse the table to get the interpolation right, for negative polarity (np.interp always @@ -124,7 +129,8 @@ def abs_first_item_or_nan(tup: Tuple): def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[complex], q_s: float, b_direct_ref: float, - b_cross_ref: float, distribution: str = 'gaussian', + b_cross_ref: float, b_direct_add: float = 0, b_cross_add: float = 0, + distribution: str = 'gaussian', fraction_of_qs_allowed_on_positive_side: float = 0.05, tolerance=1e-10): """ @@ -136,6 +142,10 @@ def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[comple the x plane or $\alpha_y \sigma_y$ if working in the y plane) :param b_cross_ref: the cross detuning coefficient multiplied by sigma (i.e. $\alpha_{xy} \sigma_y$ if working in the x plane or $\alpha_{yx} \sigma_x$ if working in the y plane) + :param b_direct_add: an additional contribution to the direct detuning coefficient which is kept fixed in the + procedure to find the stability threshold + :param b_cross_add: an additional contribution to the cross detuning coefficient which is kept fixed in the + procedure to find the stability threshold :param distribution: the transverse distribution of the beam. It can be 'gaussian' or 'parabolic' :param fraction_of_qs_allowed_on_positive_side: to determine azimuthal mode number l_mode (around which is drawn the stability diagram), one can consider positive tuneshift up to this fraction of q_s (default=5%) @@ -148,7 +158,7 @@ def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[comple # find max octupole current required from a list of modes, given their tuneshifts b_coefficients = np.array([find_detuning_coeffs_threshold( tune_shift=tune_shift, q_s=q_s, b_direct_ref=b_direct_ref, - b_cross_ref=b_cross_ref, distribution=distribution, + b_cross_ref=b_cross_ref, b_direct_add=b_direct_add, b_cross_add=b_cross_add, distribution=distribution, fraction_of_qs_allowed_on_positive_side=fraction_of_qs_allowed_on_positive_side, tolerance=tolerance) for tune_shift in tune_shifts if tune_shift is not np.nan]) From a6a2cb820a4f0b223d330edcf899ca20dd0fc7fd Mon Sep 17 00:00:00 2001 From: Lorenzo Giacomel Date: Wed, 20 Sep 2023 15:45:49 +0200 Subject: [PATCH 2/8] test: add a test for constant term in detuning coefficients --- test/test_landau_damping.py | 42 +++++++++++++++++++++++++++++++++++++ 1 file changed, 42 insertions(+) diff --git a/test/test_landau_damping.py b/test/test_landau_damping.py index cbe82b48..2a8d6485 100644 --- a/test/test_landau_damping.py +++ b/test/test_landau_damping.py @@ -60,6 +60,48 @@ def test_find_detuning_coeffs_threshold(): b_cross_ref=b_cross)[1]) +def test_find_detuning_coeffs_threshold_w_added_term(): + # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model + # test positive octupole polarity + tune_shift = -7.3622423693e-05 - 3.0188372754e-06j + q_s = 2e-3 + b_direct = 1.980315192200037e-05 + b_cross = -1.4139287608495406e-05 + + frac_add = 1/3 + + b_direct_add = b_direct * frac_add + b_cross_add = b_cross * frac_add + + assert np.isclose(b_direct * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + b_direct_ref=b_direct, + b_cross_ref=b_cross, + b_direct_add=b_direct_add, + b_cross_add=b_cross_add)[0]) + assert np.isclose(b_cross * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + b_direct_ref=b_direct, + b_cross_ref=b_cross, + b_direct_add=b_direct_add, + b_cross_add=b_cross_add)[1]) + # test negative octupole polarity + b_direct = -1.2718161244917965e-05 + b_cross = 9.08066253298482e-06 + + b_direct_add = b_direct * frac_add + b_cross_add = b_cross * frac_add + + assert np.isclose(b_direct * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + b_direct_ref=b_direct, + b_cross_ref=b_cross, + b_direct_add=b_direct_add, + b_cross_add=b_cross_add)[0]) + assert np.isclose(b_cross * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + b_direct_ref=b_direct, + b_cross_ref=b_cross, + b_direct_add=b_direct_add, + b_cross_add=b_cross_add)[1]) + + def test_find_detuning_coeffs_threshold_many_tune_shifts(): # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model tune_shift = -7.3622423693e-05 - 3.0188372754e-06j From 9ff23ae97afd2828f1347d94972687b00078e857 Mon Sep 17 00:00:00 2001 From: lgiacome Date: Tue, 26 Sep 2023 12:59:24 +0200 Subject: [PATCH 3/8] fix: use correct detuning coeffs in distance from stab diag --- pywit/landau_damping.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pywit/landau_damping.py b/pywit/landau_damping.py index f91dc0ab..a6555304 100644 --- a/pywit/landau_damping.py +++ b/pywit/landau_damping.py @@ -104,8 +104,8 @@ def f(b_direct): for t_s in b_direct_i * e * 10. ** np.arange(-3, 2, 0.01)[::e]] # note: one has to reverse the table to get the interpolation right, for negative polarity (np.interp always # wants monotonically increasing abscissae) - return tune_shift.imag - np.interp(tune_shift.real, np.real(stab)[::int(np.sign(b_direct_ref))], - np.imag(stab)[::int(np.sign(b_direct_ref))]) + return tune_shift.imag - np.interp(tune_shift.real, np.real(stab)[::int(np.sign(b_direct_i))], + np.imag(stab)[::int(np.sign(b_direct_i))]) # Newton root finding try: From 017050159ad775521920bbafec75ea429a58821e Mon Sep 17 00:00:00 2001 From: lgiacome Date: Tue, 26 Sep 2023 13:09:35 +0200 Subject: [PATCH 4/8] fix: use bisection instead of newton to find stab diagram limit --- pywit/landau_damping.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/pywit/landau_damping.py b/pywit/landau_damping.py index a6555304..c73a4d1f 100644 --- a/pywit/landau_damping.py +++ b/pywit/landau_damping.py @@ -1,6 +1,6 @@ import numpy as np from scipy.special import exp1 -from scipy.optimize import newton +from scipy.optimize import bisect from typing import Sequence, Tuple @@ -109,7 +109,11 @@ def f(b_direct): # Newton root finding try: - b_direct_new = newton(f, b_direct_ref, tol=tolerance) + # we use 1e-15 as a bound instead of 0 because 0 would cause a divsion by 0 in dispersion_integral_2d + if b_direct_ref > 0: + b_direct_new = bisect(f, 1e-15, 10 * b_direct_ref) + else: + b_direct_new = bisect(f, 100 * b_direct_ref, 1e-15) except RuntimeError: b_direct_new = np.nan else: From 73d79b829e05faf6774011c3c75aa2385137d0f6 Mon Sep 17 00:00:00 2001 From: lgiacome Date: Tue, 26 Sep 2023 13:28:13 +0200 Subject: [PATCH 5/8] style: fix several variable names --- pywit/landau_damping.py | 45 +++++++++---------- test/test_landau_damping.py | 86 ++++++++++++++++++------------------- 2 files changed, 66 insertions(+), 65 deletions(-) diff --git a/pywit/landau_damping.py b/pywit/landau_damping.py index c73a4d1f..e5ea2e03 100644 --- a/pywit/landau_damping.py +++ b/pywit/landau_damping.py @@ -59,23 +59,23 @@ def dispersion_integral_2d(tune_shift: np.ndarray, b_direct: float, b_cross: flo return -i -def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, b_direct_ref: float, b_cross_ref: float, - b_direct_add: float = 0, b_cross_add: float = 0, +def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, reference_b_direct: float, reference_b_cross: float, + added_b_direct: float = 0, added_b_cross: float = 0, fraction_of_qs_allowed_on_positive_side: float = 0.05, distribution: str = 'gaussian', tolerance=1e-10): """ Compute the detuning coefficients (multiplied by sigma) corresponding to stability diagram threshold for a complex tune shift. - It keeps fixed the ratio between b_direct_ref and b_cross_ref. + It keeps fixed the ratio between reference_b_direct and reference_b_cross. :param tune_shift: the tune shift for which the octupole threshold is computed :param q_s: the synchrotron tune - :param b_direct_ref: the direct detuning coefficient multiplied by sigma (i.e. $\alpha_x \sigma_x$ if working in + :param reference_b_direct: the direct detuning coefficient multiplied by sigma (i.e. $\alpha_x \sigma_x$ if working in the x plane or $\alpha_y \sigma_y$ if working in the y plane) - :param b_cross_ref: the cross detuning coefficient multiplied by sigma (i.e. $\alpha_{xy} \sigma_y$ if working in + :param reference_b_cross: the cross detuning coefficient multiplied by sigma (i.e. $\alpha_{xy} \sigma_y$ if working in the x plane or $\alpha_{yx} \sigma_x$ if working in the y plane) - :param b_direct_add: an additional contribution to the direct detuning coefficient which is kept fixed in the + :param added_b_direct: an additional contribution to the direct detuning coefficient which is kept fixed in the procedure to find the stability threshold - :param b_cross_add: an additional contribution to the cross detuning coefficient which is kept fixed in the + :param added_b_cross: an additional contribution to the cross detuning coefficient which is kept fixed in the procedure to find the stability threshold :param distribution: the transverse distribution of the beam. It can be 'gaussian' or 'parabolic' :param fraction_of_qs_allowed_on_positive_side: to determine azimuthal mode number l_mode (around which is drawn the @@ -93,13 +93,13 @@ def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, b_direct_ref # take away the shift from azimuthal mode number tune_shift -= q_s * l_mode - b_ratio = b_cross_ref/b_direct_ref + b_ratio = reference_b_cross / reference_b_direct if tune_shift.imag < 0.: # function to solve (distance in imag. part w.r.t stab. diagram, as a function of oct. current) def f(b_direct): - b_direct_i = b_direct + b_direct_add - b_cross_i = b_ratio * b_direct + b_cross_add + b_direct_i = b_direct + added_b_direct + b_cross_i = b_ratio * b_direct + added_b_cross stab = [-1. / dispersion_integral_2d(t_s, b_direct_i, b_cross_i, distribution=distribution) for e in (-1, 1) for t_s in b_direct_i * e * 10. ** np.arange(-3, 2, 0.01)[::e]] # note: one has to reverse the table to get the interpolation right, for negative polarity (np.interp always @@ -110,10 +110,10 @@ def f(b_direct): # Newton root finding try: # we use 1e-15 as a bound instead of 0 because 0 would cause a divsion by 0 in dispersion_integral_2d - if b_direct_ref > 0: - b_direct_new = bisect(f, 1e-15, 10 * b_direct_ref) + if reference_b_direct > 0: + b_direct_new = bisect(f, 1e-15, 10 * reference_b_direct) else: - b_direct_new = bisect(f, 100 * b_direct_ref, 1e-15) + b_direct_new = bisect(f, 100 * reference_b_direct, 1e-15) except RuntimeError: b_direct_new = np.nan else: @@ -132,23 +132,24 @@ def abs_first_item_or_nan(tup: Tuple): return np.nan -def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[complex], q_s: float, b_direct_ref: float, - b_cross_ref: float, b_direct_add: float = 0, b_cross_add: float = 0, +def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[complex], q_s: float, + reference_b_direct: float, reference_b_cross: float, + added_b_direct: float = 0, added_b_cross: float = 0, distribution: str = 'gaussian', fraction_of_qs_allowed_on_positive_side: float = 0.05, tolerance=1e-10): """ Compute the detuning coefficients corresponding to the most stringent stability diagram threshold for a sequence of - complex tune shifts. It keeps fixed the ratio between b_direct_ref and b_cross_ref. + complex tune shifts. It keeps fixed the ratio between reference_b_direct and reference_b_cross. :param tune_shifts: the sequence of complex tune shifts :param q_s: the synchrotron tune - :param b_direct_ref: the direct detuning coefficient multiplied by sigma (i.e. $\alpha_x \sigma_x$ if working in + :param reference_b_direct: the direct detuning coefficient multiplied by sigma (i.e. $\alpha_x \sigma_x$ if working in the x plane or $\alpha_y \sigma_y$ if working in the y plane) - :param b_cross_ref: the cross detuning coefficient multiplied by sigma (i.e. $\alpha_{xy} \sigma_y$ if working in + :param reference_b_cross: the cross detuning coefficient multiplied by sigma (i.e. $\alpha_{xy} \sigma_y$ if working in the x plane or $\alpha_{yx} \sigma_x$ if working in the y plane) - :param b_direct_add: an additional contribution to the direct detuning coefficient which is kept fixed in the + :param added_b_direct: an additional contribution to the direct detuning coefficient which is kept fixed in the procedure to find the stability threshold - :param b_cross_add: an additional contribution to the cross detuning coefficient which is kept fixed in the + :param added_b_cross: an additional contribution to the cross detuning coefficient which is kept fixed in the procedure to find the stability threshold :param distribution: the transverse distribution of the beam. It can be 'gaussian' or 'parabolic' :param fraction_of_qs_allowed_on_positive_side: to determine azimuthal mode number l_mode (around which is drawn the @@ -161,8 +162,8 @@ def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[comple """ # find max octupole current required from a list of modes, given their tuneshifts b_coefficients = np.array([find_detuning_coeffs_threshold( - tune_shift=tune_shift, q_s=q_s, b_direct_ref=b_direct_ref, - b_cross_ref=b_cross_ref, b_direct_add=b_direct_add, b_cross_add=b_cross_add, distribution=distribution, + tune_shift=tune_shift, q_s=q_s, reference_b_direct=reference_b_direct, + reference_b_cross=reference_b_cross, added_b_direct=added_b_direct, added_b_cross=added_b_cross, distribution=distribution, fraction_of_qs_allowed_on_positive_side=fraction_of_qs_allowed_on_positive_side, tolerance=tolerance) for tune_shift in tune_shifts if tune_shift is not np.nan]) diff --git a/test/test_landau_damping.py b/test/test_landau_damping.py index 2a8d6485..e6b812b8 100644 --- a/test/test_landau_damping.py +++ b/test/test_landau_damping.py @@ -46,18 +46,18 @@ def test_find_detuning_coeffs_threshold(): q_s = 2e-3 b_direct = 1.980315192200037e-05 b_cross = -1.4139287608495406e-05 - assert np.isclose(b_direct, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, b_direct_ref=b_direct, - b_cross_ref=b_cross)[0]) - assert np.isclose(b_cross, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, b_direct_ref=b_direct, - b_cross_ref=b_cross)[1]) + assert np.isclose(b_direct, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, + reference_b_cross=b_cross)[0]) + assert np.isclose(b_cross, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, + reference_b_cross=b_cross)[1]) # test negative octupole polarity b_direct = -1.2718161244917965e-05 b_cross = 9.08066253298482e-06 - assert np.isclose(b_direct, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, b_direct_ref=b_direct, - b_cross_ref=b_cross)[0]) - assert np.isclose(b_cross, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, b_direct_ref=b_direct, - b_cross_ref=b_cross)[1]) + assert np.isclose(b_direct, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, + reference_b_cross=b_cross)[0]) + assert np.isclose(b_cross, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, + reference_b_cross=b_cross)[1]) def test_find_detuning_coeffs_threshold_w_added_term(): @@ -68,38 +68,38 @@ def test_find_detuning_coeffs_threshold_w_added_term(): b_direct = 1.980315192200037e-05 b_cross = -1.4139287608495406e-05 - frac_add = 1/3 - - b_direct_add = b_direct * frac_add - b_cross_add = b_cross * frac_add - - assert np.isclose(b_direct * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross, - b_direct_add=b_direct_add, - b_cross_add=b_cross_add)[0]) - assert np.isclose(b_cross * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross, - b_direct_add=b_direct_add, - b_cross_add=b_cross_add)[1]) + added_fraction = 1/3 + + added_b_direct = b_direct * added_fraction + added_b_cross = b_cross * added_fraction + + assert np.isclose(b_direct * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + reference_b_direct=b_direct, + reference_b_cross=b_cross, + added_b_direct=added_b_direct, + added_b_cross=added_b_cross)[0]) + assert np.isclose(b_cross * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + reference_b_direct=b_direct, + reference_b_cross=b_cross, + added_b_direct=added_b_direct, + added_b_cross=added_b_cross)[1]) # test negative octupole polarity b_direct = -1.2718161244917965e-05 b_cross = 9.08066253298482e-06 - b_direct_add = b_direct * frac_add - b_cross_add = b_cross * frac_add + added_b_direct = b_direct * added_fraction + added_b_cross = b_cross * added_fraction - assert np.isclose(b_direct * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross, - b_direct_add=b_direct_add, - b_cross_add=b_cross_add)[0]) - assert np.isclose(b_cross * (1-frac_add), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross, - b_direct_add=b_direct_add, - b_cross_add=b_cross_add)[1]) + assert np.isclose(b_direct * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + reference_b_direct=b_direct, + reference_b_cross=b_cross, + added_b_direct=added_b_direct, + added_b_cross=added_b_cross)[0]) + assert np.isclose(b_cross * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + reference_b_direct=b_direct, + reference_b_cross=b_cross, + added_b_direct=added_b_direct, + added_b_cross=added_b_cross)[1]) def test_find_detuning_coeffs_threshold_many_tune_shifts(): @@ -112,17 +112,17 @@ def test_find_detuning_coeffs_threshold_many_tune_shifts(): b_direct = 3.960630384598084e-05 b_cross = -2.8278575218404585e-05 assert np.isclose(b_direct, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross)[0]) + reference_b_direct=b_direct, + reference_b_cross=b_cross)[0]) assert np.isclose(b_cross, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross)[1]) + reference_b_direct=b_direct, + reference_b_cross=b_cross)[1]) # test negative octupole polarity b_direct = -2.5436322491034757e-05 b_cross = 1.816132506682559e-05 assert np.isclose(b_direct, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross)[0]) + reference_b_direct=b_direct, + reference_b_cross=b_cross)[0]) assert np.isclose(b_cross, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - b_direct_ref=b_direct, - b_cross_ref=b_cross)[1]) + reference_b_direct=b_direct, + reference_b_cross=b_cross)[1]) From 671b04375c6860d7cf79897d2c3de6a476400ea3 Mon Sep 17 00:00:00 2001 From: lgiacome Date: Tue, 26 Sep 2023 14:16:46 +0200 Subject: [PATCH 6/8] test: compact tests using parametrize --- test/test_landau_damping.py | 142 ++++++++++-------------------------- 1 file changed, 40 insertions(+), 102 deletions(-) diff --git a/test/test_landau_damping.py b/test/test_landau_damping.py index e6b812b8..d38685f9 100644 --- a/test/test_landau_damping.py +++ b/test/test_landau_damping.py @@ -1,128 +1,66 @@ import pytest +from pytest import mark import numpy as np from pywit.landau_damping import dispersion_integral_2d, find_detuning_coeffs_threshold from pywit.landau_damping import find_detuning_coeffs_threshold_many_tune_shifts -def test_dispersion_integral_2d(): - - distribution = 'gaussian' - tune_shift = -7.3622423693e-05-3.0188372754e-06j +@mark.parametrize('distribution, expected_value', + [['gaussian', 7153.859519171599 - 2722.966574677259j], + ['parabolic', 6649.001778623455 - 3168.4257879737443j], + ]) +def test_dispersion_integral_2d(distribution: str, expected_value: complex): + # reference values obtained with DELPHI + tune_shift = -7.3622423693e-05 - 3.0188372754e-06j b_direct = 7.888357197059519e-05 b_cross = -5.632222163778799e-05 - # reference value obtained with DELPHI - expected_value = 7153.859519171599-2722.966574677259j - - assert np.isclose(np.real(expected_value), np.real(dispersion_integral_2d(tune_shift=tune_shift, - b_direct=b_direct, - b_cross=b_cross, - distribution=distribution))) - assert np.isclose(np.imag(expected_value), np.imag(dispersion_integral_2d(tune_shift=tune_shift, - b_direct=b_direct, - b_cross=b_cross, - distribution=distribution))) - - distribution = 'parabolic' + test_value = dispersion_integral_2d(tune_shift=tune_shift, b_direct=b_direct, b_cross=b_cross, + distribution=distribution) - # reference value obtained with DELPHI - expected_value = 6649.001778623455-3168.4257879737443j + assert np.isclose(np.real(expected_value), np.real(test_value)) + assert np.isclose(np.imag(expected_value), np.imag(test_value)) - assert np.isclose(np.real(expected_value), np.real(dispersion_integral_2d(tune_shift=tune_shift, - b_direct=b_direct, - b_cross=b_cross, - distribution=distribution))) - assert np.isclose(np.imag(expected_value), np.imag(dispersion_integral_2d(tune_shift=tune_shift, - b_direct=b_direct, - b_cross=b_cross, - distribution=distribution))) - -def test_find_detuning_coeffs_threshold(): +@mark.parametrize('b_direct, b_cross, added_b_direct, added_b_cross', + [[1.980315192200037e-05, -1.4139287608495406e-05, 0, 0], + [-1.2718161244917965e-05, 9.08066253298482e-06, 0, 0], + [1.980315192200037e-05, -1.4139287608495406e-05, 1.980315192200037e-05 / 3, + -1.4139287608495406e-05 / 3], + [-1.2718161244917965e-05, 9.08066253298482e-06, -1.2718161244917965e-05 / 3, + 9.08066253298482e-06 / 3], + ]) +def test_find_detuning_coeffs_threshold(b_direct: float, b_cross: float, added_b_direct: float, added_b_cross: float): # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model # test positive octupole polarity tune_shift = -7.3622423693e-05 - 3.0188372754e-06j q_s = 2e-3 - b_direct = 1.980315192200037e-05 - b_cross = -1.4139287608495406e-05 - assert np.isclose(b_direct, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, - reference_b_cross=b_cross)[0]) - assert np.isclose(b_cross, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, - reference_b_cross=b_cross)[1]) - # test negative octupole polarity - b_direct = -1.2718161244917965e-05 - b_cross = 9.08066253298482e-06 - assert np.isclose(b_direct, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, - reference_b_cross=b_cross)[0]) - assert np.isclose(b_cross, find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, reference_b_direct=b_direct, - reference_b_cross=b_cross)[1]) + b_direct_test, b_cross_test = find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, + reference_b_direct=b_direct, + reference_b_cross=b_cross, + added_b_direct=added_b_direct, + added_b_cross=added_b_cross) + + assert np.isclose(b_direct - added_b_direct, b_direct_test) + assert np.isclose(b_cross - added_b_cross, b_cross_test) -def test_find_detuning_coeffs_threshold_w_added_term(): +@mark.parametrize('b_direct, b_cross', + [[3.960630384598084e-05, -2.8278575218404585e-05], + [-2.5436322491034757e-05, 1.816132506682559e-05], + ]) +def test_find_detuning_coeffs_threshold_many_tune_shifts(b_direct: float, b_cross: float): # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model - # test positive octupole polarity tune_shift = -7.3622423693e-05 - 3.0188372754e-06j q_s = 2e-3 - b_direct = 1.980315192200037e-05 - b_cross = -1.4139287608495406e-05 - - added_fraction = 1/3 + tune_shifts = [np.nan, tune_shift, 2 * tune_shift] - added_b_direct = b_direct * added_fraction - added_b_cross = b_cross * added_fraction + b_direct_test, b_cross_test = find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, + reference_b_direct=b_direct, + reference_b_cross=b_cross) + assert np.isclose(b_direct, b_direct_test) + assert np.isclose(b_cross, b_cross_test) - assert np.isclose(b_direct * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross, - added_b_direct=added_b_direct, - added_b_cross=added_b_cross)[0]) - assert np.isclose(b_cross * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross, - added_b_direct=added_b_direct, - added_b_cross=added_b_cross)[1]) - # test negative octupole polarity - b_direct = -1.2718161244917965e-05 - b_cross = 9.08066253298482e-06 - added_b_direct = b_direct * added_fraction - added_b_cross = b_cross * added_fraction - - assert np.isclose(b_direct * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross, - added_b_direct=added_b_direct, - added_b_cross=added_b_cross)[0]) - assert np.isclose(b_cross * (1-added_fraction), find_detuning_coeffs_threshold(tune_shift=tune_shift, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross, - added_b_direct=added_b_direct, - added_b_cross=added_b_cross)[1]) - - -def test_find_detuning_coeffs_threshold_many_tune_shifts(): - # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model - tune_shift = -7.3622423693e-05 - 3.0188372754e-06j - q_s = 2e-3 - tune_shifts = [np.nan, tune_shift, 2*tune_shift] - - # test positive octupole polarity - b_direct = 3.960630384598084e-05 - b_cross = -2.8278575218404585e-05 - assert np.isclose(b_direct, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross)[0]) - assert np.isclose(b_cross, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross)[1]) - # test negative octupole polarity - b_direct = -2.5436322491034757e-05 - b_cross = 1.816132506682559e-05 - assert np.isclose(b_direct, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross)[0]) - assert np.isclose(b_cross, find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, - reference_b_direct=b_direct, - reference_b_cross=b_cross)[1]) From 9db91ed5d5937ffd0e465cdf8455b80dc7c5171f Mon Sep 17 00:00:00 2001 From: lgiacome Date: Tue, 26 Sep 2023 15:16:24 +0200 Subject: [PATCH 7/8] fix: wrong bounds of bisect --- pywit/landau_damping.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pywit/landau_damping.py b/pywit/landau_damping.py index e5ea2e03..1f91704a 100644 --- a/pywit/landau_damping.py +++ b/pywit/landau_damping.py @@ -111,9 +111,10 @@ def f(b_direct): try: # we use 1e-15 as a bound instead of 0 because 0 would cause a divsion by 0 in dispersion_integral_2d if reference_b_direct > 0: - b_direct_new = bisect(f, 1e-15, 10 * reference_b_direct) + b_direct_new = bisect(f, 1e-15, 100 * reference_b_direct) else: - b_direct_new = bisect(f, 100 * reference_b_direct, 1e-15) + b_direct_new = bisect(f, 100 * reference_b_direct, -1e-15) + except RuntimeError: b_direct_new = np.nan else: From c86db7c3e53416a6aad4d7d82d3ce5a440f9fb66 Mon Sep 17 00:00:00 2001 From: nimounet Date: Wed, 24 Sep 2025 16:48:23 +0200 Subject: [PATCH 8/8] Allowing both bisect and Newton root finding algorithm in threshold finding function (for Landau damping) --- pywit/landau_damping.py | 35 ++++++++++++++++++++++------------- test/test_landau_damping.py | 23 ++++++++++++++++++----- 2 files changed, 40 insertions(+), 18 deletions(-) diff --git a/pywit/landau_damping.py b/pywit/landau_damping.py index 1f91704a..c759574d 100644 --- a/pywit/landau_damping.py +++ b/pywit/landau_damping.py @@ -1,6 +1,6 @@ import numpy as np from scipy.special import exp1 -from scipy.optimize import bisect +from scipy.optimize import bisect, newton from typing import Sequence, Tuple @@ -62,7 +62,7 @@ def dispersion_integral_2d(tune_shift: np.ndarray, b_direct: float, b_cross: flo def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, reference_b_direct: float, reference_b_cross: float, added_b_direct: float = 0, added_b_cross: float = 0, fraction_of_qs_allowed_on_positive_side: float = 0.05, - distribution: str = 'gaussian', tolerance=1e-10): + distribution: str = 'gaussian', algorithm: str = 'newton', tolerance=1e-10): """ Compute the detuning coefficients (multiplied by sigma) corresponding to stability diagram threshold for a complex tune shift. @@ -80,10 +80,11 @@ def find_detuning_coeffs_threshold(tune_shift: complex, q_s: float, reference_b_ :param distribution: the transverse distribution of the beam. It can be 'gaussian' or 'parabolic' :param fraction_of_qs_allowed_on_positive_side: to determine azimuthal mode number l_mode (around which is drawn the stability diagram), one can consider positive tune shift up to this fraction of q_s (default=5%) - :param tolerance: tolerance on difference w.r.t stability diagram, for Newton's root finding + :param algorithm: root-solving algorithm ('bisect' or 'newton') + :param tolerance: tolerance on difference w.r.t stability diagram, for root finding and for the final check that the roots are actually proper roots. :return: the detuning coefficients corresponding to the stability diagram threshold if the corresponding mode is - unstable, 0 if the corresponding mode is stable or np.nan if the threshold cannot be found (failure of Newton's + unstable, 0 if the corresponding mode is stable or np.nan if the threshold cannot be found (failure of root-solving algorithm). """ # evaluate azimuthal mode number @@ -107,13 +108,20 @@ def f(b_direct): return tune_shift.imag - np.interp(tune_shift.real, np.real(stab)[::int(np.sign(b_direct_i))], np.imag(stab)[::int(np.sign(b_direct_i))]) - # Newton root finding + # root finding try: - # we use 1e-15 as a bound instead of 0 because 0 would cause a divsion by 0 in dispersion_integral_2d - if reference_b_direct > 0: - b_direct_new = bisect(f, 1e-15, 100 * reference_b_direct) + if algorithm=='bisect': + # we use 1e-15 as a bound instead of 0 because 0 would cause a divsion by 0 in dispersion_integral_2d + if reference_b_direct > 0: + b_direct_new = bisect(f, 1e-15, 100 * reference_b_direct) + else: + b_direct_new = bisect(f, 100 * reference_b_direct, -1e-15) + + elif algorithm=='newton': + b_direct_new = newton(f, reference_b_direct, tol=tolerance) + else: - b_direct_new = bisect(f, 100 * reference_b_direct, -1e-15) + raise ValueError("algorithm must be either 'newton' or 'bisect'") except RuntimeError: b_direct_new = np.nan @@ -138,7 +146,7 @@ def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[comple added_b_direct: float = 0, added_b_cross: float = 0, distribution: str = 'gaussian', fraction_of_qs_allowed_on_positive_side: float = 0.05, - tolerance=1e-10): + algorithm: str = 'newton', tolerance = 1e-10): """ Compute the detuning coefficients corresponding to the most stringent stability diagram threshold for a sequence of complex tune shifts. It keeps fixed the ratio between reference_b_direct and reference_b_cross. @@ -155,17 +163,18 @@ def find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts: Sequence[comple :param distribution: the transverse distribution of the beam. It can be 'gaussian' or 'parabolic' :param fraction_of_qs_allowed_on_positive_side: to determine azimuthal mode number l_mode (around which is drawn the stability diagram), one can consider positive tuneshift up to this fraction of q_s (default=5%) - :param tolerance: tolerance on difference w.r.t stability diagram, for Newton's root finding + :param algorithm: root-solving algorithm ('bisect' or 'newton') + :param tolerance: tolerance on difference w.r.t stability diagram, for root finding and for the final check that the roots are actually proper roots. :return: the detuning coefficients corresponding to the most stringent stability diagram threshold for all the given tune shifts if the corresponding mode is unstable, 0 if all modes are stable or np.nan if the - no threshold can be found (failure of Newton's algorithm). + no threshold can be found (failure of root-solving algorithm). """ # find max octupole current required from a list of modes, given their tuneshifts b_coefficients = np.array([find_detuning_coeffs_threshold( tune_shift=tune_shift, q_s=q_s, reference_b_direct=reference_b_direct, reference_b_cross=reference_b_cross, added_b_direct=added_b_direct, added_b_cross=added_b_cross, distribution=distribution, fraction_of_qs_allowed_on_positive_side=fraction_of_qs_allowed_on_positive_side, - tolerance=tolerance) for tune_shift in tune_shifts if tune_shift is not np.nan]) + algorithm=algorithm, tolerance=tolerance) for tune_shift in tune_shifts if tune_shift is not np.nan]) return max(b_coefficients, key=abs_first_item_or_nan) diff --git a/test/test_landau_damping.py b/test/test_landau_damping.py index d38685f9..d3feadf5 100644 --- a/test/test_landau_damping.py +++ b/test/test_landau_damping.py @@ -1,5 +1,5 @@ import pytest -from pytest import mark +from pytest import mark, raises import numpy as np from pywit.landau_damping import dispersion_integral_2d, find_detuning_coeffs_threshold @@ -23,6 +23,15 @@ def test_dispersion_integral_2d(distribution: str, expected_value: complex): assert np.isclose(np.imag(expected_value), np.imag(test_value)) +def test_wrong_algo_in_find_detuning_coeffs_threshold(): + with raises(ValueError, match="algorithm must be either 'newton' or 'bisect'"): + _ = find_detuning_coeffs_threshold(tune_shift=1e-5-1j*1e-5, q_s=2e-3, + reference_b_direct=1e-15, + reference_b_cross=1e-15, + algorithm='newto') + + +@mark.parametrize('algorithm',['newton','bisect']) @mark.parametrize('b_direct, b_cross, added_b_direct, added_b_cross', [[1.980315192200037e-05, -1.4139287608495406e-05, 0, 0], [-1.2718161244917965e-05, 9.08066253298482e-06, 0, 0], @@ -31,7 +40,8 @@ def test_dispersion_integral_2d(distribution: str, expected_value: complex): [-1.2718161244917965e-05, 9.08066253298482e-06, -1.2718161244917965e-05 / 3, 9.08066253298482e-06 / 3], ]) -def test_find_detuning_coeffs_threshold(b_direct: float, b_cross: float, added_b_direct: float, added_b_cross: float): +def test_find_detuning_coeffs_threshold(algorithm: str, b_direct: float, + b_cross: float, added_b_direct: float, added_b_cross: float): # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model # test positive octupole polarity tune_shift = -7.3622423693e-05 - 3.0188372754e-06j @@ -41,17 +51,19 @@ def test_find_detuning_coeffs_threshold(b_direct: float, b_cross: float, added_b reference_b_direct=b_direct, reference_b_cross=b_cross, added_b_direct=added_b_direct, - added_b_cross=added_b_cross) + added_b_cross=added_b_cross, + algorithm=algorithm) assert np.isclose(b_direct - added_b_direct, b_direct_test) assert np.isclose(b_cross - added_b_cross, b_cross_test) +@mark.parametrize('algorithm',['newton','bisect']) @mark.parametrize('b_direct, b_cross', [[3.960630384598084e-05, -2.8278575218404585e-05], [-2.5436322491034757e-05, 1.816132506682559e-05], ]) -def test_find_detuning_coeffs_threshold_many_tune_shifts(b_direct: float, b_cross: float): +def test_find_detuning_coeffs_threshold_many_tune_shifts(algorithm: str, b_direct: float, b_cross: float): # reference values obtained with the old impedance wake model https://gitlab.cern.ch/IRIS/HLLHC_IW_model tune_shift = -7.3622423693e-05 - 3.0188372754e-06j q_s = 2e-3 @@ -59,7 +71,8 @@ def test_find_detuning_coeffs_threshold_many_tune_shifts(b_direct: float, b_cros b_direct_test, b_cross_test = find_detuning_coeffs_threshold_many_tune_shifts(tune_shifts=tune_shifts, q_s=q_s, reference_b_direct=b_direct, - reference_b_cross=b_cross) + reference_b_cross=b_cross, + algorithm=algorithm) assert np.isclose(b_direct, b_direct_test) assert np.isclose(b_cross, b_cross_test)