Python-реализация алгоритмов и инструментарий для моделирования динамики сверхпроводникого квантового интерферометра с двумя джозефсоновскими переходами (СКВИД постоянного тока)#
1. Описание модели#
Эффект Джозефсона и джозефсоновский переход:
Связь двух сверхпроводящих слоев посредством тонкого слоя несверхпроводящего барьера образует структуру называемую джозефсоновским переходом (в честь Английского ученого Брайана Джозефсона). При пропускании электрического тока через джозефсоновский переход в зависимости от величины тока наблюдается стационарный и нестационарный эффект Джозефсона.
Стационарный эффект Джозефсона. При пропускании тока ниже критического значения (\(I<I_{c}\)) в джозефсоновском переходе (ДП) отсутсвует напряжение и через переход течет сверхпроводящий ток. Данный ток пропорционален синусу разности фаз параметров порядка сверхпроводящих слоев образующих ДП.
Нестационарный эффект Джозефсона. При увеличении тока выше критического значения (\(I>I_{c}\)) возникает переменное напряжение в ДП и оно пропорционально производной разности фаз
Сверхпроводниковые квантовые интерференционные устройства (СКВИД)
Системы содержащие джозефсоновские переходы давно нашли свое применение в различных отраслях науки и техники. Одним из таких примеров является так называемые сверхпроводниковые квантовые интерференционные устройства. В литературе эти устройства известны как СКВИД (от английского SQUID - Superconducting Quantum Interference Device). Различают двух типов СКВИДа: СКВИД постоянного тока (DC-SQUID) и СКВИД переменного тока (RF-SQUID). СКВИД постоянного тока представляет собой сверхпроводниковое кольцо, состоящее из двух параллельно соединенных джозефсоновских переходов, и широко применяется как устройство для измерения очень слабых магнитных полей. А СКВИД переменного тока применяются как чувствительный прибор слабого переменного поля.
Рассмотрим, СКВИД постоянного тока (DC-SQUID), схематический вид которого представлен на Рис.1. Плоскость на которой лежат джозефсоновские переходы и ограничено соединяющими сверхпроводящими проводниками будем называть поверхностью СКВИДа. Предположим, что перпендикулярно к поверхности СКВИДа приложено постоянное магнитное поле \(B\) направленное снизу вверх и оно приведет к потоку \(\Phi_{ext}=BS\), где \(S\) площадь поверхности СКВИДа.
Рис.1 Схематический вид СКВИДа
Рис.2 Эквивалентная схема СКВИДа.
Динамика перехода описывается Resitively Capasitevily Shunted Junction (RCSJ) моделью\cite{r}. В рамках этой модели ДП моделируется как параллельное соединение конденсатора, резистора и сверхпроводника, через которые протекают, соответственно, ток смещения, квазичастичный и сверхпроводящий токи. Учитывая индуктивности сверхпроводящих проводников соединяющие джозефсоновские переходы 1 и 2 Можно написать эквивалентную схему СКВИДа, которая представлена на Рис. 2. Индексами 1 и 2 обозначены, соответсвенно левый и правый и правый ДП. Через СКВИД пропускается внешный ток \(I\), которая распределясь на \(I_{1}\) и \(I_{2}\) течет через джозефсоновские переходы 1 и 2 и в резистивном состоянии ДП-ов приводит к среднему напряжению \(V\). Кроме того из за наличия потока внешнего магнитного поля \(\Phi_{ext}\) возникает циркулирующий сверхпроводящий ток \(J\). Здесь \(L_{1}\) и \(L_{2}\) обозначают индуктивность сверхпроводящих проводников соединяющих джозефсоновские переходы и \(C_{1,2}\), \( R_{1,2}\), \( I_{c,1,2}\) емкость, сопротивление и критический ток ДП, соответвенно. Согласно RCSJ-модели токи \(I_{1}\) и \(I_{2}\) определяются выражениями
где \(V_{i}=\frac{\Phi_{0}}{2\pi}\frac{d\varphi_{i}}{dt}\) напряжение в \(i\)-ом переходе, \(\Phi_{0}=\frac{\hbar}{2e}\)-квант магнитного потока. С другой стороны, согласно правила Кирхгофа для токов можно записать
Далее для простоты предполагаем, что индуктивности проводников равны друг другу и ДП СКВИДа одинаковы, т.е. полагаем
Магнитный поток через СКВИД складывается из потока создаваемого внешним полем \(\Phi_{ext}\) и циркулирующим током \(J\) и записывается в виде
Хорошо известно, что магнитный поток через сверхпроводниковые контуры квантуется, т.е. принимает дискретное значения кратное к \(\Phi_{0}\). Выражение для квантования магнитного потока для рассматриваемой системы можно записать
или
где \(\varphi_{ext}\) безразмерный поток внешнего магнитного поля нормированный на \(\Phi_{0}\). Отсюда можно получить выражению для циркулирующего тока
Тогда с учетом (2)-(5) можно получить замкнутую систему уравнений для описания СКВИДа:
Переходя к безразмерным величинам можно написать
где \(\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}\) - безразмерный поток внешнего магнитного поля. Здесь напряжение нормировано на \(V_{c}=I_{c}R\), ток на \(I_{c}\) и время \(t\) на \(\omega^{-1}_{c}\), где \(\omega_{c}=2\pi V_{c}/\Phi_{0}\). Начальные условия для системы уравнений \eqref{eqsystem} заданы в виде \(V_1(0)=V_2(0)=0\) и \(\varphi_1(0)=\varphi_2(0)=0\).
Постановка задачи: Необходимо вычислить:
Вольт-амперную характеристику (ВАХ) СКВИДа
Зависимость критического тока от потока внешнего магнитного поля.
Алгоритм рассчета
Ниже опишем алгоритм вычисления ВАХ и зависимости критического тока от внешнего магнитного поля. При заданном значений параметров и начальных условий при фиксированном значении тока \(I\) численно решается СДУ (\ref{eqsystem}) в интервале времени [0, Tmax] с шагом \(\Delta t\). Полученная временная зависимость напряжения усредняется по времени по формуле
и в результате вычисляется одна точка ВАХ для заданного значения \(I\). Отметим, что при усреднении интегрирования начинается от времени \(Tmin\), после которого решение стабилизируется. Далее меняется значение тока на δI и повторяется вышеописанные действия. При этом вычисленные φ(Tmax) и V (Tmax) для предыдущего значения тока будут использоваться в качестве начального условия для текущего значения тока. Таким образом увеличивая ток до Imax и уменьшая его обратно до нуля получаем ВАХ ДП.
Для определения критического тока в процессе вычисления ВАХ фиксируется значения внешнего тока при которой выполняется условие V>0, т е значение тока при которой возникает напряжение, считается критическим током.
Подключаем библиотеки#
import numpy as np
import matplotlib.pyplot as plt
import math as m
from scipy.integrate import solve_ivp
from functools import partial
from scipy.integrate import odeint
import time
from numbalsoda import lsoda_sig, lsoda, dop853
from numba import njit, cfunc
import seaborn as sns
sns.set()
sns.set(style="whitegrid")
%matplotlib inline
1. Численное решение задачи Коши для системы уравнений (1, 2)#
Функция для вычисления значений правых частей уравнений (1, 2) (Numba version)#
@cfunc(lsoda_sig)
def rhs(t, S, dS, p):
''' Функция определяющая правую часть систему уравнений.
Возвращает значения правой части уравнений при заданных значениях
параметров.
S - массив математических функций
dS - массив производных от математических функций по времени
p - массив параметров
'''
# Задаем параметров модели
beta_c = p[0]
beta_L = p[1]
Pi = p[2]
n = p[3]
Icurr = p[4]
ph_ext = p[5]
# 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*Pi*(n-ph_ext)))
dph2 = V2
dV2 = (1./beta_c) * ((Icurr/2.) - V2 - np.sin(ph2)
- (1./(2.*beta_L)) * (ph2-ph1+2*Pi*(n-ph_ext)))
# Передаем правые части уравнений в массив dS
dS[0] = dph1
dS[1] = dV1
dS[2] = dph2
dS[3] = dV2
funcptr = rhs.address # address to ODE function
Задаем значения параметров модели для численного решения системы уравнений#
Важно: необходимо соблюдать типы переменных, когда задаем начальные условия и параметры
beta_c = 0.5 # Параметр диссипации
beta_L = 1.0 # Амплитуда внешнего излучения
Pi = 3.14159265358979323
n = 0 # Период внешнего излучения
Icurr = 2
ph_ext = 0.5
Tmax = 800 # Максимальное значение времени
deltat = 0.05
t0 = 0
nt = int(Tmax / deltat)
Создаем массив с значениями времени для решения уравнения и задаем начальные условия
time_array = np.linspace(0, Tmax, num=nt)
s0 = np.array([0, 0, 0, 0]) # Initial conditions
Создаем массив с значениями параметров
# beta_c, beta_L, Pi, n, Icurr, ph_ext
data = np.array([beta_c, beta_L, Pi, n, Icurr, ph_ext])
data
array([0.5 , 1. , 3.14159265, 0. , 2. ,
0.5 ])
Численно решаем систему уравнений методом LSODA
%%time
# integrate with lsoda method
usol, success = lsoda(funcptr, s0, time_array, rtol = 1e-8, atol = 1e-8, data = data)
CPU times: user 267 ms, sys: 1.56 ms, total: 268 ms
Wall time: 268 ms
Извлекаем из массива решений значения напряжения V1 и V2
V1time = usol[:, 1]
V2time = usol[:, 3]
Временная зависимость V1 и V2
fig = plt.figure(figsize=(12, 10))
plt.plot(time_array, V1time, label='Временная зависимость', linewidth=2)
plt.xlabel('time', size=12)
plt.ylabel('V1', size=12)
plt.legend(loc='upper right')
plt.show()
fig = plt.figure(figsize=(12, 10))
plt.plot(time_array, V2time, label='Временная зависимость', linewidth=2)
plt.xlabel('time', size=12)
plt.ylabel('V2', size=12)
plt.legend(loc='upper right')
plt.show()
2. Алгоритм для нахождения вольт-амперной характеристики (ВАХ) СКВИДа#
Для нахождения ВАХ нам нужно полученную временную зависимость усреднить. Для этого необходимо усреднить по отдельности V1 и V2 и затем вычислить среднее от полученных < V1 > и < V2 >. Этот шаг можно реализовать в виде функции которая принимая временные зависимости V1 и V2 возвращает среднее значения < V >
Функция для усреднения напряжения
@njit
def averageV_new(ntmin, nt, deltat, V1, V2):
''' Функция для вычисления среднего значения напряжения.
Возвращает среднее значение напряжения для заданной
временной зависимости напряжений в первом и втором ДП
ntmin-номер шага по времени соответсвующая
нижней границе интегрирования,
nt - количество точек по времени, в которых найдены решения,
V1 - временная зависимость напряжения первого ДП,
V2 - временная зависимость напряжения второго ДП,
'''
V1av = np.sum(V1[ntmin: nt])/((nt-ntmin))
V2av = np.sum(V2[ntmin: nt])/((nt-ntmin))
Vav = (V1av + V2av)/2
return Vav
Функция для вычисления одной точки ВАХ
Создаем функцию, которая вычисляет одну точку ВАХ, т.е. выполняет следующие шаги:
принимает начальные условия и значения параметров модели
решает численно уравнения
используя ранее созданную функцию averageV_new усредняет напряжение
возвращает среднее напряжение и начальные условия для вычисления следующей точки ВАХ
@njit
def cvcpoint_new(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat):
''' Функция для вычисления одной точки ВАХ.
Возвращает для заданного значения тока вычисленную среднее значение
напряжения при заданных значений параметров.
s0 - массив начальных условий,
beta_c, beta_L, Pi, n, Icurr, ph_ext,-параметры модели
nt - количество точек по времени, в которых найдены решения,
ntmin - номер шага по времени соответсвующее
нижней границе интегрирования,
t0 - начальное значение времени,
deltat - шаг по времени
'''
data = np.array([beta_c, beta_L, Pi, n, Icurr, ph_ext])
t_e = np.linspace(t0, nt*deltat, nt)
# integrate with lsoda method
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]
s0 = np.array([ph1time[nt-1], V1time[nt-1], ph2time[nt-1], V2time[nt-1]])
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, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI):
''' Функция для вычисления однопетлевой ВАХ.
Возвращает для заданного значения параметров массива
с значениями тока и вычисленного напряжения
s0 - массив начальных условий,
beta_c, beta_L, Pi, n, Icurr, ph_ext,-параметры модели
nt - количество точек по времени, в которых найдены решения,
ntmin - номер шага по времени соответсвующее
нижней границе интегрирования,
t0 - начальное значение времени,
deltat - шаг по времени
deltaI - шаг по току
'''
Vplot = []
Iplot = []
while Icurr < 100:
res = cvcpoint_new(s0, beta_c, beta_L, Pi, 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)))
Задаем значения параметров для вычисления ВАХ
Tmin = 50.0 # Начало интервала для интегрирования для усреднения
ntmin = int(Tmin/deltat)
deltaI = 0.001
Icurr = 0.0
a = 1
Imax = 2.2
s0 = np.array([0.0, 0.0, 0.0, 0.0])
res = cvcpoint_new(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat)
print(res)
print(res[0])
print(res[1:])
Вычислим ВАХ
start = time.time()
res = comput_vax(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
График ВАХ
fig = plt.figure(figsize=(12, 10))
plt.plot(res[:, 0], res[:, 1], label='ВАХ', linewidth=3.0)
plt.xlabel('I', size=12)
plt.ylabel('V', size=12)
plt.legend(loc='upper left')
plt.show()
Описание ВАХ#
Как видно из ВАХ с увеличением тока начиная от нуля до 0.8 значения напряжение равна нулю (стационарный эффект Джозефсона). При достижении значения тока 0.8 появляется напряжение и далее с ростом тока растет и напряжение. Как было отмечено в описании модели ток при котором появляется напряжение назывется критическим током. Теперь рассмотрим как критический ток меняется с изменением внешнего потока.
Исследование влияния \(\varphi_{ext}\) на ВАХ#
Проведем вычисления при различных значений \(\varphi_{ext}\), например при \(\varphi_{ext}=0\), \(\varphi_{ext}=0.25\) и \(\varphi_{ext}=0.5\). Численное моделирование проводилось в интервале времени [0, 1000] с шагом 0.05.
ph_ext = 0
start = time.time()
res0 = comput_vax(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
ph_ext = 0.25
start = time.time()
res1 = comput_vax(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
ph_ext = 0.5
start = time.time()
res2 = comput_vax(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat, a, deltaI)
end = time.time()
print(end-start)
График ВАХ#
fig = plt.figure(figsize=(9, 8))
plt.plot(res0[:, 0], res0[:, 1], label='ВАХ phi_ext=0', linewidth=3.0)
plt.plot(res1[:, 0], res1[:, 1], label='ВАХ phi_ext=0.25', linewidth=3.0)
plt.plot(res2[:, 0], res2[:, 1], label='ВАХ phi_ext=0.5', linewidth=3.0)
plt.xlabel('I', size=12)
plt.ylabel('V', size=12)
plt.legend(loc='upper left')
plt.savefig('cvcsquid.png')
plt.show()
Как видно из графика ВАХ при \(\varphi_{ext}=0\) критический ток равен \(I_{c}=2\), при \(\varphi_{ext}=0.25\) - \(I_{c}=1.58\), а при \(\varphi_{ext}=0.5\) - \(I_{c}=0.8\), т.е. с изменением \(\varphi_{ext}\) меняется величина критического тока.
3. Вычисление зависимости критического тока от внешнего магнитного поля#
Создаем фукнцию getIc.#
Для этого вводим следующий параметр
epsilon - точность критического тока
ph_ext = 0.5
epsilon = 0.001
Icurr = 0
@njit
def get_Ic(s0, beta_c, beta_L, Pi, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon):
''' Функция для вычисления значении критического тока
при заданных параметров. Возвращает для заданного
значения параметров значение критического тока.
s0 - массив начальных условий,
beta_c, beta_L, Pi, n, Icurr, ph_ext,-параметры модели
nt - количество точек по времени, в которых найдены решение,
ntmin - номер шага по времени соответсвующее
нижней границе интегрирования,
t0-начальное значение времени,
deltat - шаг по времени
deltaI - шаг по току
epsilon - точность нахождения критического тока
'''
Icurr = 0
a = 1
while Icurr < 100:
res = cvcpoint_new(s0, beta_c, beta_L, Pi, n,
Icurr, ph_ext, nt, ntmin, t0, deltat)
Vav = res[0]
s0 = res[1:]
if(Vav > epsilon):
Ic = Icurr
break
Icurr += deltaI
return Ic
epsilon = 0.01
ph_ext = 1
t_start = time.time()
k = get_Ic(s0, beta_c, beta_L, Pi, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon)
t_finish = time.time()
print(f'Execution time {t_finish - t_start} s')
print(k)
Исследование зависимости критического тока от внешнего магнитного поля#
Для этого вводим следующие параметры:
npoint - количество точек в интервале значений внешнего поля
ph_extmin - минимальное значение внешнего поля
ph_extmax - максимальное значение внешнего поля
deltaph_ext - шаг изменения внешнего поля
npoint = 51
ph_extmin = 0.0
ph_extmax = 2.0
deltaph_ext = (ph_extmax - ph_extmin) / (npoint-1)
print(deltaph_ext)
Создаем массив значений внешнего поля для которых нужно вычислить критическй ток
ph_ext_array = np.linspace(ph_extmin, ph_extmax, num=npoint)
ph_ext_array.shape
(51,)
Создаем пустой массив для значений крического тока
Ic_array = np.zeros(npoint)
Ic_array.shape
(51,)
Создаем функцию зависимости критического тока от внешнего поля#
@njit
def funk_numba(j, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon):
''' Функция для вычисления значении критического тока
при заданных параметров и фиксированном значении
потока внешнего магнитного поля с номером j
в интервале [ph_extmin, j*deltaph_ext]. Возвращает для
заданного значения ph_ext значение критического тока.
s0 - массив начальных условий,
beta_c, beta_L, Pi, n, Icurr, ph_ext,-параметры модели
nt - количество точек по времени, в которых найдены решение,
ntmin - номер шага по времени соответсвующее
нижней границе интегрирования,
t0-начальное значение времени,
deltat - шаг по времени
deltaI - шаг по току
epsilon - точность критического тока
'''
ph_ext = ph_extmin + deltaph_ext * j
Ic = get_Ic(s0, beta_c, beta_L, Pi, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon)
return Ic
epsilon = 0.1
%%time
for i in range(0, npoint):
Icrit = funk_numba(i, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n, ph_ext,
nt, ntmin, t0, deltat, deltaI, epsilon)
Ic_array[i] = Icrit
print(Ic_array)
Посторим график зависимости критического тока от внешнего поля#
fig = plt.figure(figsize=(6, 6))
plt.plot(ph_ext_array, Ic_array,
label='Зависимость $I_c$ от phi_ext', linewidth=3.0)
plt.xlabel('$phi_{ext}$', size=12)
plt.ylabel('$I_{c}$', size=12)
plt.legend(loc='upper right')
plt.savefig('Icfield.png')
plt.show()
Данный график демонстрирует модулирование критического тока с изменением \(\varphi_{ext}\).
Параллелизация с 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'
1:80: E501 line too long (81 > 79 characters)
config.NUMBA_DEFAULT_NUM_THREADS
88
set_num_threads(5)
config.NUMBA_NUM_THREADS = 20
config.NUMBA_NUM_THREADS
20
config.THREADING_LAYER = 'omp'
threading_layer()
'tbb'
config.THREADING_LAYER
'omp'
# demonstrate the threading layer chosen
print("Threading layer chosen: %s" % threading_layer())
get_num_threads()
5
set_num_threads(5)
get_num_threads()
5
import os
print('NUMBA NUM THREADS')
try:
print(os.environ['NUMBA_NUM_THREADS'])
except:
print(config.NUMBA_DEFAULT_NUM_THREADS)
5:1: E722 do not use bare 'except'
import numba as nb
from numba import prange
@njit(parallel=True)
def comput_parallel(npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n,
ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon):
''' Функция для вычисления зависимости критического тока
от потока внешнего магнитного поля при заданных параметров.
Возвращает массив значений критических токов для заданного
интервала значений ph_ext.
npoint - количество точек ph_ext в интервале
ph_extmin - начало интервала потока внешнего магнитного поля,
deltaph_ext - шаг потока внешнего магнитного поля,
s0 - массив начальных условий,
beta_c, beta_L, Pi, n, Icurr, ph_ext,-параметры модели
nt - количество точек по времени, в которых найдены решение,
ntmin - номер шага по времени соответсвующее
нижней границе интегрирования,
t0-начальное значение времени,
deltat - шаг по времени
deltaI - шаг по току
epsilon - точность критического тока
'''
Ic_array = np.empty(npoint)
Ic_array[:] = 0
for i in nb.prange(npoint):
Ic = funk_numba(i, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n,
ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon)
Ic_array[i] = Ic
return Ic_array
%%time
Ic_array = comput_parallel(npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n, ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon)
fig = plt.figure(figsize=(6, 6))
plt.plot(ph_ext_array, Ic_array,
label='Зависимость Ic от phi_ext', linewidth=3.0)
plt.xlabel('$phi_{ext}$', size=12)
plt.ylabel('$I_{c}$', size=12)
plt.legend(loc='upper right')
plt.show()
Вычисления при разных чисел потоков.
На этом этапе реализована параллелизация вычислении критического тока. На графике представлена рассчитанная зависимость критического тока СКВИДа от потока \(\varphi_{ext}\) в интервале потока [0, 2]. Данная зависимость демонстрирует модулирование критического тока с изменением \(\varphi_{ext}\), что согласуется с известными результатами. Теперь можно вычислить зависимость времени вычисления от числа вычислительных потоков.
set_num_threads(1)
t_start = time.time()
Ic_array = comput_parallel(npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n, ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon)
t_finish = time.time()
print(f'Execution time {t_finish - t_start} s')
Execution time 51.1058554649353 s
set_num_threads(5)
t_start = time.time()
Ic_array = comput_parallel(npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n, ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon)
t_finish = time.time()
print(f'Execution time {t_finish - t_start} s')
Execution time 13.315101623535156 s
set_num_threads(10)
t_start = time.time()
Ic_array = comput_parallel(npoint, ph_extmin, deltaph_ext, s0, beta_c, beta_L, Pi, n, ph_ext, nt, ntmin, t0, deltat, deltaI, epsilon)
t_finish = time.time()
print(f'Execution time {t_finish - t_start} s')
Execution time 7.867552280426025 s
Построим график зависимости времени вычисления от числа потоков
Proc=[1,5,10,15,20,25,30,35,40,45,50,55,60,65,70]
Time=[51.1,13.31,7.86,5.65,4.9,5.05,8.5,8.53,9.9,9.81,11.91,11.35,11.37,11.65,11.56]
fig = plt.figure(figsize=(6, 6))
plt.plot(Proc, Time, label='Calculation time from number of Threads', linewidth=3.0)
plt.xlabel('Threads number', size=12)
plt.ylabel('Time (s)', size=12)
plt.legend(loc='upper right')
plt.savefig('calc_time.png')
plt.show()
Разработанный инструментарий позволяет моделировать физические свойства сверхпроводникового квантового интерферометра с двумя джозефсоновскими переходами (СКВИД постоянного тока). Реализованы алгоритмы для вычисления вольт-амперной характеристики СКВИДа под воздействием внешнего магнитного поля и зависимости критического тока СКВИДа от внешнего магнитного поля. При программной реализации разработанных алгоритмов использовались функции библиотеки Numba, в том числе механизм распараллеливания при вычислении зависимости критического тока от потока внешнего магнитного поля. Представленные результаты были получены на экосистеме ML/DL/HPC Гетерогенной платформы HybriLIT (ЛИТ ОИЯИ).
Работа выполнена при финансовой поддержки РНФ в рамках проекта No 22-71-10022.