Python-реализация алгоритмов и инструментарий для моделирования динамики сверхпроводникого квантового интерферометра с двумя джозефсоновскими переходами (СКВИД постоянного тока)#
Блокнот предназначен моделировани физических свойств сверхпроводящего квантового интерференционного устройства (СКВИД) с двумя джозефсоновскими пеереходами (СКВИД постоянного тока) и исследования зависимости величины тока возврата от величины потока внешнего магнитного поля, а также демонстрации возникновения модуляции напряжения от величины внешнего магнитного поля.
Разработаны алгоритмы для вычисления вольт-амперной характеристики СКВИДа под действием внешнего магнитного поля, зависимосней тока возврата СКВИД от внешнего магнитного поля и демонстрации возникновения модуляции напряжения СКВИДа от величины внешнего магнитного поля при постоянном токе. В программной реализации разработанных алгоритмов использовались функции библиотеки Numba, в том числе и механизм распараллеливания вычислений между вычислительными ядрами многоядерных CPU при вычислении зависимостей тока возврата от внешнего поля и модуляции напряжения от потока внешнего поля при постоянном токе.
Представленные результаты были получены на экосистеме ML/DL/HPC Гетерогенной платформы HybriLIT (ЛИТ ОИЯИ).
Работа выполнена при финансовой поддержки РНФ в рамках проекта No 22-71-10022.
1. Математическая модель СКВИДа#
Эффект Джозефсона и джозефсоновский переход
Связь двух сверхпроводящих слоев, посредством тонкого слоя несверхпроводящего барьера образует структуру, называемую джозефсоновским переходом (в честь английского ученого Брайана Джозефсона). При пропускании электрического тока через джозефсоновский переход, в зависимости от величины тока, наблюдаются стационарный и нестационарный эффекты Джозефсона [1],[2].
Стационарный эффект Джозефсона: при пропускании внешнего тока ниже критического значения (\(I < I_{c}\)), в джозефсоновском переходе (ДП) отсутсвует напряжение и через ДП протрекает сверхпроводящий ток. Данный ток пропорционален синусу разности фаз параметров порядка сверхпроводящих слоев образующих ДП:
Нестационарный эффект Джозефсона: при пропускании внешнего тока выше критического значения (\(I > I_{c}\)), возникает переменное напряжение в ДП, и оно пропорционально производной разности фаз:
Согласно \(RCSJ\) модели МакКамбера-Стюарта [3],[4], ДП-ы можно представить в виде параллельных контактов резистора, сверхпроводника и конденсатора, как показано на Рис. 1., а ток, протекающий ДП, состоит из суммы квазичастичного тока \(\frac{V}{R}\), сверхпроводящего тока \(I_{c}\sin\varphi\) и тока смещения \(C\frac{dV}{dt}\). Здесь \(R\) и \(C\) - соответственно, сопротивление и емкость ДП в резистивном состоянии.
Рис. 1. Эквивалентная схема джозефсоновского перехода. |
Сверхпроводниковые квантовые интерференционные устройства (СКВИД)
Рассмортрим сверхпроводниковое квантовое интерференционное устройство СКВИД (SQUID) [5],[6], состоящее из двух параллельно соединенных сверхпроводящим кольцом джозефсоновских переходов - двуконтактный СКВИД []. Предполагается, что через СКВИД пропускается электрический ток и сверхпроводящее кольцо СКВИДа находится во внешнем магнитном поле, как показано на Рис. 2., а эквивалентная схема СКВИД представлена на Рис. 3.
Рис. 2. Схема двуконтактного СКВИДа. |
Рис. 3. Эквивалентная схема двуконтактного СКВИДа. |
В частном случае предположим, что индуктивности проводников равны(\(L_{1}=L_{2}=L\)), а также джозефсоновские переходы одинаковы, т.е. их сопротивления, емкости и значения их критических токов соответственно равны (\(R_{1}=R_{2}=R, \quad C_{1}=C_{2}=C, \quad I_{c1}=I_{c2}=I_{c}\)).
Полный ток, протекающий через СКИД складывается из суммы токов, протекающих через каждый ДП, и циркулирующего сверхпроводящего тока \(J\), вызванного потоком внешнего магнитного поля.
Циркулирующий сверхпроводящий ток задается выражением:
где \(\Phi_{0}\) - квант магнитного потока, \(B\)- внешнее магнитное поле, пронизывающее кольцо СКВИДа, \(S\) - эффективная площадь, ограниченная сверхпроводящим кольцом и двумя джозефсоновскими переходами, \(n\) - целое число, \(\varphi_{1}\) и \(\varphi_{2}\) - разности фаз параметров порядка сверхпроводников, соответственно, первого и второго ДП-ов.
Согласно правилу Кирхгофа, токи, протекающие через каждый ДП, можно записать следующими выражениями:
где \(I\) - полный ток СКВИДа.
Математическая модель СКВИД-а, состоящего из двух одинаковых ДП, находящегося во внешнем магнитном поле, записывается согласно правиду Кирхгофа и \(RCSJ\) модели джозефсоновских переходов:
Переходя к безразмерным величинам, система уравнений (5) записывается в виде:
Здесь напряжение нормировано на \(V_{c}=I_{c}R\), ток нормирован на \(I_{c}\) и время \(t\) нормировано на \(\omega^{-1}_{c}\), а
\(\beta_{c}=C\Phi_{0}\omega_{c}^{2}/(2\pi I_{c})\) - параметр МакКамбера,
\(\omega_{c}=2\pi R I_{c}/\Phi_{0}\)- характеристическая частота Джозефсона,
\(\beta_{L}=2\pi L I_{c}/\Phi_{0}\)- индуктивность в безразмерных единицах,
\(\Phi_{0}=h/(2e)\) - квант магнитного потока,
\(\varphi_{ext}=\Phi_{ext}/\Phi_{0}\) - безразмерный поток внешнего магнитного поля,
\(\omega_{c}=2\pi V_{c}/\Phi_{0}\) - характеристическая частота джозефсона.
Таким образом, система уравнений (6) определяет задачу Коши, для которой начальные условия задаются в виде: \(V_1(0)=V_2(0)=0\) и \(\varphi_1(0)=\varphi_2(0)=0\).
Влияние параметра МакКамбера \(\beta_c\)
Ключевым безразмерным параметром, определяющим режим работы джозефсоновского перехода [7], является параметр МакКамбера \(\beta_c = \omega_c R C\).
В зависимости от величины \(\beta_c\) возможны два качественно различных режима:
Режим сильного затухания (overdamped) при \(\beta_c \ll 1\). В этом случае емкостные эффекты пренебрежимо малы (\(C \to 0\)), и динамика перехода описывается моделью RSJ (Resistively Shunted Junction). Вольт-амперная характеристика (ВАХ) такого перехода не имеет гистерезиса: переход в резистивное состояние и обратно происходит при одном и том же значении тока, то есть ток возврата \(I_r\) практически равен критическому току \(I_c\) (см. Рис. 3а).
Режим слабого затухания (underdamped) при \(\beta_c \gg 1\). Здесь емкость играет существенную роль, накапливая заряд и создавая инерционность. Вольт-амперная характеристика такого джозефсоновского перехода обладает ярко выраженным гистерезисом: переключение в резистивное состояние происходит при критическом токе \(I_c\), а обратное переключение в сверхпроводящее состояние — при значительно меньшем токе возврата \(I_r < I_c\) (см. Рис. 3б). Именно это расхождение между \(I_c\) и \(I_r\) является отличительной чертой таких джозефсоновских переходов.
Рис. 3a. ВАХ СКВИД при $\beta_{C}=0.01$
Рис. 3б. ВАХ СКВИД при $\beta_{C}=25$Таким образом, варьируя параметр \(\beta_c\) (например, изменяя емкость перехода \(C\) или его сопротивление \(R\)), можно получать различные виды ВАХ как с гистерезисом \(\left(\Delta I = I_c - I_r \right)\) и без него.
Постановка задачи#
Необходимо рассчитать:
Вольт-амперные характеристики (ВАХ) при различных значениях внешнего магнитного потока \(\varphi_{ext}\) для случаев \(\beta_c \ll 1\) и \(\beta_c \gg 1\).
Зависимость величины критического тока (для \(\beta_c \ll 1\)) и тока возврата (для \(\beta_c \gg 1\)) от внешнего магнитного потока \(\varphi_{ext}\).
Зависимость модуляции напряжения на СКВИДе от внешнего магнитного потока \(\varphi_{ext}\) в режиме сильного затухания (\(\beta_c \ll 1\)).
Алгоритм нахождения ВАХ и величины тока возврата
Ниже опишем алгоритм нахождения ВАХ.
При заданном значений параметров модели СКВИДа и начальных условиях численно решается задача Коши для системы уравнений [6] в интервале времени \([0, T_{max}]\) с шагом \(\Delta t\) при фиксированном значении тока \(I\).
Полученные для каждого ДП временные зависимости напряжений усредняется по времени по формуле:
\[ \bar{V}_{i} = \frac{1}{T_{\max} - T_{\min}} \int_{T_{\min}}^{T_{\max}} V_{i}(t) \, dt, \quad i=1,2. \tag{7} \]Среднее напряжение на СКВИД вычисляется по формуле:
\[ \bar{V} = \frac{\bar{V}_{1} + \bar{V}_{2}}{2}. \tag{8} \]Таким образом вычисляется одна точка ВАХ для заданного значения \(I\). Отметим, что при усреднении, интегрирование начинается от времени \(T_{min}\), которое выбирается таким образом, чтобы исключить влияние переходных процессов, т.е. решение выходит на стационарный режим. При этом для обеспечения непрерывности ВАХ (особенно в гистерезисной области) в качестве начальных условий для нового значения тока используются конечные значения \(\varphi_{i}(T_{\max})\) и \(V_{i}(T_{\max})\), полученные при предыдущем значении тока.
Последовательно увеличивая ток от \(I=0\) до \(I_{\max}\) (прямой ход), а затем уменьшая его обратно до нуля (обратный ход), получаем полную ВАХ СКВИДа, включая возможную гистерезисную ветвь.
Для определения тока возврата \(I_r\) в процессе вычисления ВАХ на ветви уменьшения тока фиксируется такое значение внешнего тока \(I\), при котором усредненное напряжение на СКВИДе \(\bar{V}\) становится равным нулю (с заданной точностью). Это значение и есть искомая величина тока возврата \(I_r\) для данного внешнего магнитного поля \(\varphi_{ext}\) при \(\beta_c \gg 1\).
Список литературы#
1. B.D. Josephson, Possible new effects in superconductive tunnelling, J. Phys. Lett. 1, 251-253 (1962). DOI: 10.1016/0031-9163(62)91369-0;
2. Anderson P. W., Rowell J. M., Probable Observation of the Josephson Superconducting Tunneling Effect, Physical Review Letters. 1963. Vol. 10, iss. 6. P. 230-232. DOI: 10.1103/PhysRevLett.10.230;
3. D. E. McCumber J., Effect of ac Impedance on dc Voltage‐Current Characteristics of Superconductor Weak‐Link Junctions, Appl. Phys. 39, 3114 (1968). DOI: 10.1063/1.1656743;
4. Stewart W. C., Current-Voltage Characteristics of Josephson Junctions, Applied Physics Letters. 1968. Vol. 12, iss. 8. P. 277-280. DOI: 10.1063/1.1651991;
5. R. C. Jaklevic, J. Lambe, A. H. Silver, J. E. Mercereau «Quantum Interference from a Static Vector Potential in a Superconducting Field», Phys. Rev. Lett. 12, 159 (1964). DOI: 10.1103/PhysRevLett.12.159;
6. Clarke J., Braginski A.I. (Eds.) The SQUID Handbook. Vol. I. Fundamentals and Technology of SQUIDs and SQUID Systems. - Weinheim: Wiley-VCH, 2004. - 408 p. - ISBN 3-527-40229-2;
7. R. Grimaudo, D. Valenti, B. Spagnolo, A. Troisi, G. Filatrella, C. Guarcello «Axion field influence on Josephson junction quasipotential»,Phys. Rev. B 108, 174521 (2023). DOI: 10.1103/PhysRevB.108.174521.
2. Численное решение задачи Коши для системы уравнений (6)#
Для численного решения поставленных задач и моделирования динамики СКВИД добавим следующие библиотеки и модули:
NumPy и math — для выполнения векторных и математических операций;
Matplotlib и seaborn — для визуализации результатов;
numbalsoda (lsoda, dop853) — для численного интегрирования систем обыкновенных дифференциальных уравнений (ОДУ);
Numba (njit, cfunc) — для ускорения вычислений путём JIT-компиляции правых частей уравнений;
time — для оценки производительности алгоритмов.
Matplotlib - для построения графических зависимостей,
seaborn — для расширенной визуализации, что позволилит повысить наглядность представления результатов.
Данный набор инструментов обеспечивает как высокую точность, так и приемлемое время рассчёта.
import numpy as np
import matplotlib.pyplot as plt
import math as m
from functools import partial
import time
from numbalsoda import lsoda_sig, lsoda, dop853
from numba import njit, cfunc
import seaborn as sns
sns.set()
sns.set(style="whitegrid", context="talk", font_scale=1.2)
%matplotlib inline
2.1. Алгоритм вычисления временных зависимостей напряжения на джозефсоновских переходах#
Задаем функцию для вычисления значений правых частей уравнений (6) (Numba version)
@cfunc(lsoda_sig)
def rhs(t, S, dS, p):
''' Функция определяющая правую часть системы уравнений.
Возвращает значения правой части уравнений при заданных параметрах.
S - массив математических функций;
dS - массив производных от математических функций по времени;
p - массив параметров СКВИДа;
t - временной параметр.
'''
# Задаем значения параметров модели
beta_c = p[0]
beta_L = p[1]
n = p[2]
Icurr = p[3]
ph_ext = p[4]
# S массив, который содержит значения функций
ph1 = S[0]
V1 = S[1]
ph2 = S[2]
V2 = S[3]
# Задаем правую часть системы уравнений
dph1 = V1
dV1 = (1./beta_c) * ((Icurr/2.) - V1 - np.sin(ph1)
+ (1./(2.*beta_L)) * (ph2-ph1+2*np.pi*(n-ph_ext)))
dph2 = V2
dV2 = (1./beta_c) * ((Icurr/2.) - V2 - np.sin(ph2)
- (1./(2.*beta_L)) * (ph2-ph1+2*np.pi*(n-ph_ext)))
# Передаем правые части уравнений в массив dS
dS[0] = dph1
dS[1] = dV1
dS[2] = dph2
dS[3] = dV2
funcptr = rhs.address # указатель ODE функции
Создаем массив значений параметров модели
beta_c = 0.1 # Параметр МакКамбера
beta_L = 1.0 # Безразмерная индуктивность
n = 0.0 # Число квантов потока
Icurr = 2.0 # Значение внешнего тока I
ph_ext = 0.5 # Значение потока внешнего магнитного поля
Создаем массив data, в который записываем параметры модели
data = np.array([beta_c, beta_L, n, Icurr, ph_ext])
data
array([0.1, 1. , 0. , 2. , 0.5])
Создаем равномерную сетку по времени \(\left[0, T_{max} \right]\) в виде массива
t0 = 0 # Начальное значение времени
Tmax = 150 # Максимальное значение времени
deltat = 0.05 # Величина шага по времени
nt = int((Tmax-t0) / deltat) # Количечство значений в интервале времени
time_array = np.linspace(0, Tmax, num=nt)
Создаем массив начальных условий
Задаем значения \(\varphi_1(t_{0})\), \(V_1(t_{0})\), \(\varphi_2(t_{0})\), \(V_2(t_{0})\) в начальный момент времени
s0 = np.array([0, 0, 0, 0])
Численно решаем задачу Коши методом LSODA
%%time
usol, success = lsoda(funcptr, s0, time_array, rtol = 1e-8, atol = 1e-8, data = data)
CPU times: user 109 ms, sys: 940 μs, total: 110 ms
Wall time: 110 ms
Извлекаем из массива решений значения напряжения \(V_1(t)\) и \(V_2(t)\)
V1time = usol[:, 1]
V2time = usol[:, 3]
Строим графики зависимости напряжения на каждом ДП от времени:
fig, axes = plt.subplots(2, 1, figsize=(8, 8))
axes[0].plot(time_array, V1time, 'b-', linewidth=1.75)
axes[0].set_xlabel('time')
axes[0].set_ylabel('$V_1$')
axes[0].grid(True, alpha=0.3)
axes[1].plot(time_array, V2time, 'r-', linewidth=1.75)
axes[1].set_xlabel('time')
axes[1].set_ylabel('$V_2$')
axes[1].grid(True, alpha=0.3)
fig.suptitle('Рис. 4. Временные зависимости напряжений\n'
'на джозефсоновских переходах',
fontsize=16, y=0.02)
plt.tight_layout()
plt.show()
2.2. Алгоритм нахождения вольт-амперной характеристики (ВАХ) СКВИДа#
Для вычисления ВАХ необходимо временные зависимости напряжения \(V_1(t)\) и \(V_2(t)\) усреднить по отдельности и вычислить среднее от полученных \(\bar{V_1}\) и \(\bar{V_2}\). В результате получим среднее значение напряжение \(\bar{V}\) СКВИД-а.
Функция усреднения напряжения
@njit
def averageV_new(ntmin, nt, deltat, V1, V2):
''' Функция для вычисления среднего значения напряжения.
Возвращает среднее значение напряжения для заданных
временных зависимостей напряжения на первом и втором ДП;
ntmin-номер шага по времени, соответсвующий
нижнему пределу интегрирования;
nt - количество точек по времени, в которых найдены решения;
V1 - временная зависимость напряжения первого ДП;
V2 - временная зависимость напряжения второго ДП.
'''
Vav1 = np.mean(V1[ntmin: nt])
Vav2 = np.mean(V2[ntmin: nt])
Vav = (Vav1+Vav2)/2
return Vav
Создаем функцию, которая вычисляет одну точку ВАХ, т.е. выполняет следующие операции:
принимает начальные условия и значения параметров модели;
решает численно уравнения систему уравнений (6);
используя ранее созданную функцию, averageV_new, усредняет напряжение;
возвращает среднее напряжение и начальные условия для вычисления следующей точки ВАХ.
@njit
def cvcpoint_new(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat):
''' Функция для вычисления одной точки ВАХ.
Возвращает для заданного значения тока вычисленное среднее значение
напряжения при заданных значениях параметров.
s0 - массив начальных условий;
beta_c, beta_L, n, Icurr, ph_ext - параметры модели;
nt - количество точек по времени, в которых найдены решения;
ntmin - номер шага по времени, соответсвующий
нижнему пределу интегрирования;
t0 - начальное значение времени;
deltat - шаг по времени.
'''
data = np.array([beta_c, beta_L, n, Icurr, ph_ext])
t_e = np.linspace(t0, nt*deltat, nt)
# численное интегрирование методом lsoda
usol, success = lsoda(funcptr, s0, t_e, rtol=1e-8, atol=1e-8, data=data)
# извлечение результатов из массива usol
ph1time = usol[:, 0]
V1time = usol[:, 1]
ph2time = usol[:, 2]
V2time = usol[:, 3]
# вызов функции нахождения среднего напряжения
Vav = averageV_new(ntmin, nt, deltat, V1time, V2time)
return np.array([Vav, ph1time[nt-1], V1time[nt-1],
ph2time[nt-1], V2time[nt-1]])
Функция вычисления ВАХ:
@njit
def comput_vax(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI):
''' Функция для вычисления однопетлевой ВАХ.
Возвращает для заданного значения параметров массив
значений тока и вычисленного напряжения.
s0 - массив начальных условий;
beta_c, beta_L, n, Icurr, ph_ext - параметры модели;
nt - количество точек в интервале времени, в которых найдены решения;
ntmin - номер шага по времени, соответсвующий
нижнему пределу интегрирования;
t0 - начальное значение времени;
deltat - шаг по времени;
deltaI - шаг по току.
'''
Vplot = []
Iplot = []
while Icurr < 100:
res = cvcpoint_new(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat)
Vav = res[0]
s0 = res[1:]
Vplot.append(Vav)
Iplot.append(Icurr)
Icurr += a*deltaI
if (Icurr > Imax):
a = -1
if ((Icurr < 0) and (a == -1)):
break
return np.column_stack((np.array(Iplot), np.array(Vplot)))
Задаем значения параметров модели СКВИД для вычисления ВАХ
beta_c = 25.0 # Параметр МакКамбера
beta_L = 1.0 # Номированная индуктивность
n = 0 # Число квантов потока
ph_ext = 0.0 # Значение потока внешнего магнитного поля
Tmin = 50.0 # Начало интервала для интегрирования для усреднения
ntmin = int(Tmin/deltat) # Порядковый номер значения времени начала усреднения
deltaI = 0.0005 # Величина шага в интервале по току
Icurr = 0.0 # Начальное значение тока
Imax = 2.7 # Максимальное значение тока
a = 1 # Параметр направления расчета
s0 = np.array([0.0, 0.0, 0.0, 0.0]) # Массив начальных условий
Вычисляем ВАХ
start = time.time()
res = comput_vax(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
3.832197427749634
Строим графики ВАХ для случаев увеличения (прямой ход) и уменьшения (обратный ход) тока.
%matplotlib inline
fig, axes = plt.subplots(2, 1, figsize=(8, 8))
axes[0].plot(np.array_split(res[:, 0], 2)[0], np.array_split(res[:, 1], 2)[0],
'b-', label='$V(I_{up})$', linewidth=3.0)
axes[0].set_ylabel('$V$')
axes[0].set_xlabel('I')
axes[0].grid(True, alpha=0.5)
axes[0].legend(loc='upper left', fontsize=12)
axes[0].arrow(x=1.0, y=0.1, dx=0.5, dy=0.0,
head_width=0.1, head_length=0.2,
fc='black', ec='black', linewidth=2)
axes[0].set_title('Прямой ход', fontsize=14)
axes[1].plot(np.array_split(res[:, 0], 2)[1], np.array_split(res[:, 1], 2)[1],
'r-', label='$V(I_{down})$', linewidth=3.0)
axes[1].set_xlabel('I')
axes[1].set_ylabel('$V$')
axes[1].grid(True, alpha=0.5)
axes[1].legend(loc='upper left', fontsize=12)
axes[1].set_title('Обратный ход', fontsize=14)
axes[1].arrow(x=1.5, y=0.1, dx=-0.5, dy=0.0,
head_width=0.1, head_length=0.2,
fc='black', ec='black', linewidth=2)
fig.suptitle('Рис. 5. Вольт-амперная характеристика СКВИД',
fontsize=16, y=0.02)
plt.subplots_adjust(hspace=0.3)
plt.tight_layout()
plt.show()
3. Исследование влияния величины внешнего поля \(\varphi_{ext}\) на вид ВАХ#
Проведем расчеты для построения графиков вольт-амперных характеристик СКВИД при различных значениях внешнего поля: \(\varphi_{ext}=0\), \(\varphi_{ext}=0.25\), \(\varphi_{ext}=0.5\).
ph_ext = 0
start = time.time()
res0 = comput_vax(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
3.51435923576355
ph_ext = 0.25
start = time.time()
res1 = comput_vax(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
5.077096223831177
ph_ext = 0.5
start = time.time()
res2 = comput_vax(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
3.111905813217163
ВАХ при различных \(\varphi_{ext}\)
fig, axes = plt.subplots(2, 1, figsize=(8, 8), sharex=True)
axes[0].plot(np.array_split(res0[:, 0], 2)[0],
np.array_split(res0[:, 1], 2)[0], 'b-',
label='$V(I_{up})$, $\\varphi_{ext}=0$', linewidth=3.0)
axes[0].plot(np.array_split(res1[:, 0], 2)[0],
np.array_split(res1[:, 1], 2)[0], 'r-',
label='$V(I_{up})$, $\\varphi_{ext}=0.25$', linewidth=3.0)
axes[0].plot(np.array_split(res2[:, 0], 2)[0],
np.array_split(res2[:, 1], 2)[0], 'g-',
label='$V(I_{up})$, $\\varphi_{ext}=0.5$', linewidth=3.0)
for x_pos in np.arange(0.2, Imax, 0.2):
axes[0].axvline(x=x_pos, color='black', linestyle='dashed',
alpha=0.3, linewidth=1)
axes[0].set_ylabel('$V$')
axes[0].grid(True, alpha=0.5)
axes[0].legend(loc='upper left', fontsize=12)
axes[0].set_title('Прямой ход', fontsize=12)
axes[1].plot(np.array_split(res0[:, 0], 2)[1],
np.array_split(res0[:, 1], 2)[1], 'b-',
label='$V(I_{down})$, $\\varphi_{ext}=0$', linewidth=3.0)
axes[1].plot(np.array_split(res1[:, 0], 2)[1],
np.array_split(res1[:, 1], 2)[1], 'r-',
label='$V(I_{down})$, $\\varphi_{ext}=0.25$', linewidth=3.0)
axes[1].plot(np.array_split(res2[:, 0], 2)[1],
np.array_split(res2[:, 1], 2)[1], 'g-',
label='$V(I_{down})$, $\\varphi_{ext}=0.5$', linewidth=3.0)
for x_pos in np.arange(0.2, Imax, 0.2):
axes[1].axvline(x=x_pos, color='black', linestyle='dashed',
alpha=0.3, linewidth=1)
axes[1].set_xlabel('I')
axes[1].set_ylabel('$V$')
axes[1].grid(True, alpha=0.5)
axes[1].legend(loc='upper left', fontsize=12)
axes[1].set_title('Обратный ход', fontsize=12)
fig.suptitle(
'Рис. 6. Вольт-амперные характеристики СКВИД\n'
'при разных значениях $\\varphi_{ext}$',
fontsize=16, y=0.02)
plt.subplots_adjust(hspace=0.3)
plt.tight_layout()
plt.show()
Эти графики показывают, что, при изменении величины \(\varphi_{ext}\), происходит смещение ВАХ. При этом, значения критического тока и тока возврата, рассчитанные для разных значений внешнего поля \(\varphi_{ext}\), отличаются.
4. Вычисление зависимости величины тока возврата от внешнего магнитного поля#
Расчет зависимости тока возврата от внешнего магнитного поля включает расчёт вольт-амперных характеристик СКВИДа в диапазоне значений внешнего магнитного поля и последующее определение на каждой ВАХ значения тока, при котором напряжение становится равным нулю при движении от больших токов к нулю (при обратном ходе).
Для этого вводим следующие параметры:
npoint - количество точек в интервале значений внешнего поля
ph_extmin - минимальное значение внешнего поля
ph_extmax - максимальное значение внешнего поля
deltaph_ext - шаг изменения внешнего поля
npoint = 1000
ph_extmin = 0.0
ph_extmax = 2.0
deltaph_ext = (ph_extmax - ph_extmin) / (npoint-1)
Создаем равномерную сетку значений внешнего поля [ph_extmin, ph_extmax] с шагом deltaph_ext в виде массива
ph_ext_array = np.linspace(ph_extmin, ph_extmax, num=npoint)
ph_ext_array.shape
(1000,)
Создаем пустой массив, в который будем записывать полученные значения тока возврата.
Ireturn_array = np.zeros(npoint)
Ireturn_array.shape
(1000,)
Обозначим величину тока, до которой будет проводиться расчет ВАХ.
Imax = 2.2
Задаем функцию для определения значения величины тока возврата для заданных значений параметров модели
@njit
def get_Ireturn(s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon):
''' Функция для определения значения тока возврата
при заданных параметрах модели.
Возвращает значение тока возврата.
s0 - массив начальных условий,
beta_c, beta_L, n, Icurr, ph_ext - параметры модели,
nt - количество точек по времени,
в которых найдены решение,
ntmin - номер шага по времени, соответсвующий
нижнему пределу интегрирования,
t0 - начальное значение времени,
deltat - шаг по времени
deltaI - шаг по току
epsilon - точность нахождения тока возврата
'''
Icurr = Imax
a = 1
while Icurr < 100:
res = cvcpoint_new(s0, beta_c, beta_L, n,
Icurr, ph_ext, nt, ntmin, t0, deltat)
Vav = res[0]
s0 = res[1:]
Icurr += a*deltaI
if (Icurr > Imax):
a = -1
if (Vav <= epsilon) and (a == -1):
Ireturn = Icurr
break
return Ireturn
Задаем значения параметров модели и начальные условия:
beta_c = 25.0 # Параметр МакКамбера
beta_L = 1.0 # Номированная индуктивность
n = 0 # Число квантов потока
ph_ext = 0.0 # Значение потока внешнего магнитного поля
s0 = np.array([0.0, 0.0, 0.0, 0.0])
Также вводим параметр точности \(\varepsilon\), сравнением с которым будем определять величину тока возврата.
epsilon = 0.001
Для каждого значения величины внешнего магнитного поля в интервале вычисляем значение тока возврата
%%time
for i in range(0, npoint):
ph_ext = ph_extmin + deltaph_ext * i
Ireturn = get_Ireturn(s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon)
Ireturn_array[i] = Ireturn
IOStream.flush timed out
IOStream.flush timed out
CPU times: user 16min 57s, sys: 1.19 s, total: 16min 58s
Wall time: 17min 1s
print(Ireturn_array.shape, ph_ext_array.shape)
(1000,) (1000,)
График зависимости тока возврата от величины внешнего магнитного поля
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(ph_ext_array, Ireturn_array,
'b-', label='$I_{r}(\\varphi_{ext})$', linewidth=3.0)
ax.set_xlabel('$\\varphi_{ext}$')
ax.set_ylabel('$I_{r}$')
ax.grid(True, alpha=0.5)
fig.suptitle(
'Рис. 7. Зависимость величины тока возврата от внешнего магнитного поля',
fontsize=12, y=0.02)
plt.tight_layout()
plt.show()
4. Модуляция напряжения при фиксированном токе в зависимости от внешнего магнитного поля#
Для демонстрации модуляции напряжение СКВИД в зависимости от величины внешнего поля необходимо расчитать вольт-амперные характеристики СКВИДа в диапазоне значений внешнего магнитного поля и определение на каждой ВАХ значения напряжения, которые соответствуют фиксированным значениям тока.
Для этого вводим следующие параметры:
ph_ext_min - минимальное значение внешнего поля;
ph_extmax - максимальное значение внешнего поля;
n_point - количество точек в интервале внешнего поля;
delta_ph_ext - шаг изменения внешнего поля.
n_point = 300
ph_extmin = 0.0
ph_extmax = 2
delta_ph_ext = (ph_extmax - ph_extmin) / (n_point-1)
Создаем массив значений внешнего поля:
ph_ext_array = np.linspace(ph_extmin, ph_extmax, num=n_point)
ph_ext_array.shape
(300,)
Также вводим параметры:
n_fix_point - кол-во точек в интервале по току;
I_fix_min - минимальное значение тока в интервале;
I_fix_max - максимальное значение тока в интервале;
delta_I_fix - шаг в интервале по току.
n_fix_point = 20
I_fix_min = 0.2
delta_I_fix = 0.1
I_fix_max = I_fix_min + n_fix_point * delta_I_fix
Создаем массив фиксируемых значений тока:
I_fix_array = np.linspace(I_fix_min, I_fix_max, n_fix_point)
I_fix_array.shape
(20,)
Создаем фукцию, которая проводит расчет ВАХ при заданном значении внешнего магнитного поля и нахождит те значения напряжения, которые соответствуют значениям из массива фиксируемых значений тока I_fix_array.
@njit
def numba_V_mod(j, ph_extmin, deltaph_ext, s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, I_fix_array, n_fix_point):
'''Функция для расчета ВАХ и извлечения значений напряжения,
которые соответствуют фиксированным значениям тока
при заданных параметрах.
Возвращает для заданных значений параметров,
массива фиксированных токов и значений внешнего поля
массив значений напряжения, извлеченный из ВАХ.
s0 - массив начальных условий;
beta_c, beta_L, n, Icurr, ph_ext - параметры модели;
nt - количество точек по времени, в которых найдены решение;
ntmin - номер шага по времени, соответсвующий
нижнему пределу интегрирования;
t0 - начальное значение времени;
deltat - шаг по времени;
deltaI - шаг по току.
'''
V_I = []
a = 1
ph_ext = ph_extmin + j * deltaph_ext
cvc = comput_vax(s0, beta_c, beta_L, n, Icurr, ph_ext,
nt, ntmin, t0, deltat, a, deltaI)
cvc_0 = cvc[:, 0]
for i in I_fix_array:
# Поиск ближайщего значений тока ВАХ к значениям из I_fix_array
min_index = np.argmin(np.abs(cvc_0 - i))
V_I.append(cvc[min_index][1])
return V_I
Задаем значения параметров модели и начальные условия:
beta_c = 0.1 # Параметр МакКамбера
beta_L = 1.0 # Номированная индуктивность
n = 0 # Число квантов потока
Icurr = 0.0 # Начальное значение внешнего тока
ph_ext = 0.0 # Значение потока внешнего магнитного поля
s0 = np.array([0.0, 0.0, 0.0, 0.0])
Создаем массив, в который будем добавлять вычисленнные значения напряжения, соответствующие фиксируемым значениям тока для каждого значения внешнего поля
Vmod_array = np.empty((n_point, n_fix_point), dtype=np.float64)
Для каждого значения внешнего поля вычисляем ВАХ, находим те значения напряжения, которые соответствуют значениям тока из I_fix_array, и записываем их.
%%time
for i in range(0, n_point):
Vmod = numba_V_mod(i, ph_extmin, delta_ph_ext, s0, beta_c, beta_L, n,
ph_ext, nt, ntmin, t0, deltat,
deltaI, I_fix_array, n_fix_point)
Vmod_array[i] = Vmod
CPU times: user 18min 36s, sys: 1.62 s, total: 18min 38s
Wall time: 18min 41s
Строим график модуляции напряжения от внешнего магнитного поля
fig, ax = plt.subplots(figsize=(8, 8))
for a in range(0, n_fix_point):
ax.plot(ph_ext_array, Vmod_array[:, a], linewidth=3.0)
ax.set_xlabel('$\\varphi_{ext}$')
ax.set_ylabel('$V$')
fig.suptitle(
'Рис. 8. Модуляция напряжения\n'
'от величины внешнего магнитного поля',
fontsize=14, y=0.02)
ax.grid(True, alpha=0.5)
plt.tight_layout()
plt.show()
5. Параллелизация расчета с numba#
В этом пункте представлены ранее реализованные алгоритмы расчета зависимости величины тока возврата от величины внешнего поля и расчетов для демонстрации модуляции напрядения при фиксированных значениях тока в зависимости от внешнего поля.
Материалы документации Numba
Setting the Number of Threads
The number of threads used by numba is based on the number of CPU cores available (see numba.config.NUMBA_DEFAULT_NUM_THREADS), but it can be overridden with the NUMBA_NUM_THREADS environment variable.
https://numba.readthedocs.io/en/stable/developer/threading_implementation.html
from numba import config, njit, threading_layer
from numba import set_num_threads, get_num_threads
# set the threading layer before any parallel target compilation
config.THREADING_LAYER = 'threadsafe'
config.NUMBA_DEFAULT_NUM_THREADS
256
set_num_threads(5)
config.NUMBA_NUM_THREADS = 20
config.NUMBA_NUM_THREADS
20
config.THREADING_LAYER = 'omp'
threading_layer()
'omp'
config.THREADING_LAYER
'omp'
# demonstrate the threading layer chosen
print("Threading layer chosen: %s" % threading_layer())
Threading layer chosen: omp
get_num_threads()
5
set_num_threads(80)
get_num_threads()
80
5.1. Параллелизация алгоритма расчета тока возврата от внешнего поля#
import os
print('NUMBA NUM THREADS')
try:
print(os.environ['NUMBA_NUM_THREADS'])
except KeyError:
print(config.NUMBA_DEFAULT_NUM_THREADS)
NUMBA NUM THREADS
256
import numba as nb
from numba import prange
Задаем значения параметров модели и массив начальных условий
beta_c = 25.0 # Параметр МакКамбера
beta_L = 1.0 # Номированная индуктивность
n = 0 # Число квантов потока
Icurr = 0.0 # Начальное значение внешнего тока
ph_ext = 0.0 # Значение потока внешнего магнитного поля
s0 = np.array([0.0, 0.0, 0.0, 0.0])
npoint = 400
ph_extmin = 0.0
ph_extmax = 2.0
deltaph_ext = (ph_extmax - ph_extmin) / (npoint - 1)
ph_ext_array = np.linspace(ph_extmin, ph_extmax, npoint)
ph_ext_array.shape
(400,)
@njit(parallel=True)
def comput_parallel_I_return(npoint, ph_extmin, deltaph_ext, s0,
beta_c, beta_L, n, ph_ext, nt, ntmin,
t0, deltat, deltaI, epsilon):
''' Функция для определения значения тока возврата
при заданных параметрах модели.
Возвращает значение тока возврата.
s0 - массив начальных условий,
beta_c, beta_L, n, Icurr, ph_ext - параметры модели,
nt - количество точек по времени,
в которых найдены решение,
ntmin - номер шага по времени, соответсвующий
нижнему пределу интегрирования,
t0 - начальное значение времени,
deltat - шаг по времени
deltaI - шаг по току
epsilon - точность нахождения тока возврата
'''
Ireturn_array = np.zeros(npoint)
for i in nb.prange(npoint):
ph_ext = ph_extmin + i * deltaph_ext
Ireturn = get_Ireturn(s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon)
Ireturn_array[i] = Ireturn
return Ireturn_array
%%time
Ireturn_array = comput_parallel_I_return(
npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L,
n, ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon)
CPU times: user 9min 41s, sys: 108 ms, total: 9min 41s
Wall time: 9.98 s
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(ph_ext_array, Ireturn_array,
'b-', label='$I_{r}(\\varphi_{ext})$', linewidth=3.0)
ax.set_xlabel('$\\varphi_{ext}$')
ax.set_ylabel('$I_{r}$')
ax.grid(True, alpha=0.5)
fig.suptitle(
'Рис. 9. Зависимость величины тока возврата\n'
'от внешнего магнитного поля',
fontsize=12, y=0.02)
plt.tight_layout()
plt.show()
5.2. Параллелизация нахождения модуляции напряжения в зависимости от внешнего поля#
set_num_threads(80)
ph_ext = 0.0
Icurr = 2.7
ph_ext_min - минимальное значение внешнего поля;
ph_extmax - максимальное значение внешнего поля;
n_point - количество точек в интервале внешнего поля;
delta_ph_ext - шаг изменения внешнего поля.
n_point = 200
ph_extmin = 0.0
ph_extmax = 2
delta_ph_ext = (ph_extmax - ph_extmin) / (n_point)
n_fix_point - кол-во точек в интервале по току;
I_fix_min - минимальное значение тока в интервале;
I_fix_max - максимальное значение тока в интервале;
delta_I_fix - шаг в интервале по току.
n_fix_point = 100
I_fix_min = 0.2
delta_I_fix = 0.1
I_fix_max = I_fix_min + n_fix_point * delta_I_fix
I_fix_array - массив точек по току. ph_ext_array - массив точек по внешнему полю.
beta_c = 0.001 # Параметр МакКамбера
beta_L = 1.0 # Номированная индуктивность
n = 0 # Число квантов потока
Icurr = 0.0 # Начальное значение внешнего тока
ph_ext = 0.0 # Значение потока внешнего магнитного поля
s0 = np.array([0.0, 0.0, 0.0, 0.0])
I_fix_array = np.linspace(I_fix_min, I_fix_max, n_fix_point)
I_fix_array.shape
(100,)
ph_ext_array = np.linspace(ph_extmin, ph_extmax, n_point)
ph_ext_array.shape
(200,)
@njit(parallel=True)
def comput_parallel_V_mod(n_point, ph_extmin, delta_ph_ext, s0, beta_c,
beta_L, n, ph_ext, nt, ntmin, t0, deltat,
deltaI, I_fix_array, n_fix_point):
'''Функция для расчета ВАХ и извлечения значений напряжения,
которые соответствуют фиксированным значениям тока
при заданных параметрах.
Возвращает для заданных значений параметров,
массива фиксированных токов и значений внешнего поля
массив значений напряжения, извлеченный из ВАХ.
s0 - массив начальных условий;
beta_c, beta_L, n, Icurr, ph_ext - параметры модели;
nt - количество точек по времени, в которых найдены решение;
ntmin - номер шага по времени, соответсвующий
нижнему пределу интегрирования;
t0 - начальное значение времени;
deltat - шаг по времени;
deltaI - шаг по току.
'''
Vmod_array = np.empty((n_point, n_fix_point), dtype=np.float64)
for i in nb.prange(n_point):
Vmod = numba_V_mod(i, ph_extmin, delta_ph_ext, s0, beta_c,
beta_L, n, ph_ext, nt, ntmin, t0, deltat,
deltaI, I_fix_array, n_fix_point)
Vmod_array[i] = Vmod
return Vmod_array
Создаем массив, в который будем добавлять вычисленнные значения напряжения, соответствующие фиксируемым значениям тока для каждой ВАХ.
Vmod_array = []
%%time
Vmod_array = comput_parallel_V_mod(
n_point, ph_extmin, delta_ph_ext, s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, I_fix_array, n_fix_point)
CPU times: user 1h 49min 34s, sys: 482 ms, total: 1h 49min 35s
Wall time: 2min 32s
Обозначим np.array(Vmod_array) = data и вычислим значение num - количество кривых модуляции
data = np.array(Vmod_array)
num = len(data[0])
fig, ax = plt.subplots(figsize=(8, 8))
for a in range(0, n_fix_point):
ax.plot(ph_ext_array, Vmod_array[:, a], linewidth=3.0)
ax.set_xlabel('$\\varphi_{ext}$')
ax.set_ylabel('$V$')
fig.suptitle(
'Рис. 10. Модуляция напряжения\n'
'от величины внешнего магнитного поля',
fontsize=14, y=0.02)
ax.grid(True, alpha=0.5)
plt.tight_layout()
plt.show()
Вычисления для разного числа потоков.#
На этом этапе проводится сравнение получаемого ускорения вычислений от числа задействуемых в вычислениях количество параллельных потоков.
Задается массив значений количества потоков, а также массивы, в которых будет записывать время вычислений зависимости тока возврата и модуляции напряжения от внешнего поля.
Proc = [1, 5, 10, 15, 20, 25, 30,
35, 40, 45, 50, 55, 60,
65, 70, 75, 80]
Time_Ir = []
Time_Vmod = []
Проводим вычисления и для каждого числа потоков фиксируем время, затраченное на вычисление зависимости тока возврата от внешнего поля.
for a in Proc:
set_num_threads(a)
t_start = time.time()
I_return = comput_parallel_I_return(
npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon)
t_finish = time.time()
print(f'Execution time {t_finish - t_start} s')
Time_Ir.append(t_finish - t_start)
Execution time 1448.29776096344 s
Execution time 457.6115849018097 s
Execution time 285.73713874816895 s
Execution time 197.42999744415283 s
Execution time 145.70413398742676 s
Execution time 123.77684378623962 s
Execution time 109.3470983505249 s
Execution time 97.43126010894775 s
Execution time 81.43774056434631 s
Execution time 74.11320471763611 s
Execution time 66.94642424583435 s
Execution time 67.67852354049683 s
Execution time 58.605732679367065 s
Execution time 50.83831310272217 s
Execution time 50.988717794418335 s
Execution time 50.74732542037964 s
Execution time 42.52572154998779 s
Рассчитываем степень ускорения вычислений от числа потоков, сравнивая каждое время со временем вычисления в однопоточном режиме.
speedup_Ir = []
for i in Time_Ir:
a = Time_Ir[0] / i
speedup_Ir.append(a)
Строим график зависимости ускорения вычислений от числа задействуемых вычислительных потоков.
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(Proc, speedup_Ir,
label='Calculation time from number of Threads',
linewidth=4.0)
ax.set_xlabel('Threads number', size=14)
ax.set_ylabel('Speedup', size=14)
ax.grid(True, alpha=0.5)
fig.suptitle(
'Рис. 11. Ускорение вычислений\n'
'для разного числа потоков',
fontsize=14, y=0.02)
plt.show()
Аналогично рассчитывается зависимость ускорения вычислений от числа потоков для задачи вычисления модуляции напряжения от внешнего поля и строим график.
for a in Proc:
set_num_threads(a)
t_start = time.time()
Vmod_array = comput_parallel_V_mod(
n_point, ph_extmin, delta_ph_ext, s0, beta_c, beta_L, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, I_fix_array, n_fix_point)
t_finish = time.time()
print(f'Execution time {t_finish - t_start} s')
Time_Vmod.append(t_finish - t_start)
Execution time 5779.449589967728 s
Execution time 1416.7193474769592 s
Execution time 847.7308938503265 s
Execution time 587.268212556839 s
Execution time 465.3166027069092 s
Execution time 377.9883382320404 s
Execution time 333.5157036781311 s
Execution time 288.80485582351685 s
Execution time 243.06809306144714 s
Execution time 242.85298800468445 s
Execution time 196.40459775924683 s
Execution time 195.27308702468872 s
Execution time 203.5623426437378 s
Execution time 150.9199070930481 s
Execution time 152.85308980941772 s
Execution time 152.25207018852234 s
Execution time 151.6734857559204 s
speedup_Vmod = []
for i in Time_Vmod:
a = Time_Vmod[0] / i
speedup_Vmod.append(a)
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(Proc, speedup_Vmod,
label='Calculation time from number of Threads',
linewidth=4.0)
ax.set_xlabel('Threads number', size=14)
ax.set_ylabel('Speedup', size=14)
ax.grid(True, alpha=0.5)
ax.set_title('Рис. 12. Ускорение вычислений\n'
'для разного числа потоков',
fontsize=14)
plt.show()
