Python-реализация алгоритмов и инструментарий для моделирования динамики сверхпроводникого квантового интерферометра с двумя джозефсоновскими переходами (СКВИД постоянного тока)

Python-реализация алгоритмов и инструментарий для моделирования динамики сверхпроводникого квантового интерферометра с двумя джозефсоновскими переходами (СКВИД постоянного тока)#

1. Описание модели#

Эффект Джозефсона и джозефсоновский переход:

Связь двух сверхпроводящих слоев посредством тонкого слоя несверхпроводящего барьера образует структуру называемую джозефсоновским переходом (в честь Английского ученого Брайана Джозефсона). При пропускании электрического тока через джозефсоновский переход в зависимости от величины тока наблюдается стационарный и нестационарный эффект Джозефсона.

Стационарный эффект Джозефсона. При пропускании тока ниже критического значения ( $I < I_{c}$ ) в джозефсоновском переходе (ДП) отсутсвует напряжение и через переход течет сверхпроводящий ток. Данный ток пропорционален синусу разности фаз параметров порядка сверхпроводящих слоев образующих ДП.

\begin{array}{r} (1) & I_{s} = I_{c} \sin φ . \end{array}

Нестационарный эффект Джозефсона. При увеличении тока выше критического значения ( $I > I_{c}$ ) возникает переменное напряжение в ДП и оно пропорционально производной разности фаз

\begin{array}{r} (2) & V = \frac{ℏ}{2 e} \frac{d φ}{d t} . \end{array}

Сверхпроводниковые квантовые интерференционные устройства (СКВИД)

Системы содержащие джозефсоновские переходы давно нашли свое применение в различных отраслях науки и техники. Одним из таких примеров является так называемые сверхпроводниковые квантовые интерференционные устройства. В литературе эти устройства известны как СКВИД (от английского SQUID - Superconducting Quantum Interference Device). Различают двух типов СКВИДа: СКВИД постоянного тока (DC-SQUID) и СКВИД переменного тока (RF-SQUID). СКВИД постоянного тока представляет собой сверхпроводниковое кольцо, состоящее из двух параллельно соединенных джозефсоновских переходов, и широко применяется как устройство для измерения очень слабых магнитных полей. А СКВИД переменного тока применяются как чувствительный прибор слабого переменного поля.

Рассмотрим, СКВИД постоянного тока (DC-SQUID), схематический вид которого представлен на Рис.1. Плоскость на которой лежат джозефсоновские переходы и ограничено соединяющими сверхпроводящими проводниками будем называть поверхностью СКВИДа. Предположим, что перпендикулярно к поверхности СКВИДа приложено постоянное магнитное поле $B$ направленное снизу вверх и оно приведет к потоку $Φ_{e x t} = B S$ , где $S$ площадь поверхности СКВИДа.

Рис.1 Схематический вид СКВИДа

Рис.2 Эквивалентная схема СКВИДа.

Динамика перехода описывается Resitively Capasitevily Shunted Junction (RCSJ) моделью\cite{r}. В рамках этой модели ДП моделируется как параллельное соединение конденсатора, резистора и сверхпроводника, через которые протекают, соответственно, ток смещения, квазичастичный и сверхпроводящий токи. Учитывая индуктивности сверхпроводящих проводников соединяющие джозефсоновские переходы 1 и 2 Можно написать эквивалентную схему СКВИДа, которая представлена на Рис. 2. Индексами 1 и 2 обозначены, соответсвенно левый и правый и правый ДП. Через СКВИД пропускается внешный ток $I$ , которая распределясь на $I_{1}$ и $I_{2}$ течет через джозефсоновские переходы 1 и 2 и в резистивном состоянии ДП-ов приводит к среднему напряжению $V$ . Кроме того из за наличия потока внешнего магнитного поля $Φ_{e x t}$ возникает циркулирующий сверхпроводящий ток $J$ . Здесь $L_{1}$ и $L_{2}$ обозначают индуктивность сверхпроводящих проводников соединяющих джозефсоновские переходы и $C_{1, 2}$ , $R_{1, 2}$ , $I_{c, 1, 2}$ емкость, сопротивление и критический ток ДП, соответвенно. Согласно RCSJ-модели токи $I_{1}$ и $I_{2}$ определяются выражениями

\begin{array}{r} {\begin{cases} I_{1} = C \frac{d V_{1}}{d t} + \frac{V_{1}}{R_{1}} + I_{c 1} \sin φ_{1}, \\ I_{2} = C \frac{d V_{2}}{d t} + \frac{V_{2}}{R_{2}} + I_{c 2} \sin φ_{2}, \end{cases} \end{array}

где $V_{i} = \frac{Φ_{0}}{2 π} \frac{d φ_{i}}{d t}$ напряжение в $i$ -ом переходе, $Φ_{0} = \frac{ℏ}{2 e}$ -квант магнитного потока. С другой стороны, согласно правила Кирхгофа для токов можно записать

\begin{array}{r} {\begin{cases} I_{1} = \frac{I}{2} + J, \\ I_{2} = \frac{I}{2} - J . \end{cases} \end{array}

Далее для простоты предполагаем, что индуктивности проводников равны друг другу и ДП СКВИДа одинаковы, т.е. полагаем

C_{1} = C_{2} = C, R_{1} = R_{2} = R, L_{1} = L_{2} = L, I_{c 1} = I_{c 2} = I_{c} .

Магнитный поток через СКВИД складывается из потока создаваемого внешним полем $Φ_{e x t}$ и циркулирующим током $J$ и записывается в виде

Φ = Φ_{e x t} + L (I_{1} - I_{2}) = Φ_{e x t} + 2 L J .

Хорошо известно, что магнитный поток через сверхпроводниковые контуры квантуется, т.е. принимает дискретное значения кратное к $Φ_{0}$ . Выражение для квантования магнитного потока для рассматриваемой системы можно записать

φ_{2} - φ_{1} + 2 π n = \frac{2 π}{Φ_{0}} Φ = \frac{2 π}{Φ_{0}} (Φ_{e x t} + 2 L J),

или

φ_{2} - φ_{1} + 2 π n = 2 π φ_{e x t} + \frac{4 π L}{Φ_{0}} J,

где $φ_{e x t}$ безразмерный поток внешнего магнитного поля нормированный на $Φ_{0}$ . Отсюда можно получить выражению для циркулирующего тока

J = \frac{Φ_{0}}{4 π L} [φ_{2} - φ_{1} + 2 π (n - φ_{e x t})] .

Тогда с учетом (2)-(5) можно получить замкнутую систему уравнений для описания СКВИДа:

\begin{array}{r} {\begin{cases} C \frac{d V_{1}}{d t} + \frac{V_{1}}{R} + I_{c} \sin φ_{1} = \frac{I}{2} + \frac{Φ_{0}}{4 π L} [φ_{2} - φ_{1} + 2 π (n - φ_{e x t})], \\ \frac{Φ_{0}}{2 π} \frac{d φ_{1}}{d t} = V_{1} \\ C \frac{d V_{2}}{d t} + \frac{V_{2}}{R} + I_{c} \sin φ_{2} = \frac{I}{2} - \frac{Φ_{0}}{4 π L} [φ_{2} - φ_{1} + 2 π (n - φ_{e x t})], \\ \frac{Φ_{0}}{2 π} \frac{d φ_{2}}{d t} = V_{2} \end{cases}, \end{array}

Переходя к безразмерным величинам можно написать

\begin{array}{r} {\begin{cases} \frac{d V_{1}}{d t} = \frac{1}{β_{c}} {\frac{I}{2} - \frac{d φ_{1}}{d t} - \sin (φ_{1}) + \frac{1}{2 β_{L}} [φ_{2} - φ_{1} + 2 π (n - φ_{e x t})]} \\ \frac{d φ_{1}}{d t} = V_{1} \\ \frac{d V_{2}}{d t} = \frac{1}{β_{c}} {\frac{I}{2} - \frac{d φ_{2}}{d t} - \sin (φ_{2}) - \frac{1}{2 β_{L}} [φ_{2} - φ_{1} + 2 π (n - φ_{e x t})]} \\ \frac{d φ_{2}}{d t} = V_{2} \end{cases}, \end{array}

где $β_{c} = C Φ_{0} ω_{c}^{2} / (2 π I_{c})$ - параметр МакКамбера, $ω_{c} = 2 π R I_{c} / Φ_{0}$ - характеристическая частота Джозефсона, $β_{L} = 2 π L I_{c} / Φ_{0}$ - обезразмеренная индуктивность, $Φ_{0} = h / (2 e)$ - квант магнитного потока, $φ_{e x t} = Φ_{e x t} / Φ_{0}$ - безразмерный поток внешнего магнитного поля. Здесь напряжение нормировано на $V_{c} = I_{c} R$ , ток на $I_{c}$ и время $t$ на $ω_{c}^{- 1}$ , где $ω_{c} = 2 π V_{c} / Φ_{0}$ . Начальные условия для системы уравнений \eqref{eqsystem} заданы в виде $V_{1} (0) = V_{2} (0) = 0$ и $φ_{1} (0) = φ_{2} (0) = 0$ .

Постановка задачи: Необходимо вычислить:

Вольт-амперную характеристику (ВАХ) СКВИДа
Зависимость критического тока от потока внешнего магнитного поля.

Алгоритм рассчета

Ниже опишем алгоритм вычисления ВАХ и зависимости критического тока от внешнего магнитного поля. При заданном значений параметров и начальных условий при фиксированном значении тока $I$ численно решается СДУ (\ref{eqsystem}) в интервале времени [0, Tmax] с шагом $Δ t$ . Полученная временная зависимость напряжения усредняется по времени по формуле

\bar{V} = \frac{1}{T_{max} - T_{min}} \int_{T_{min}}^{T_{max}} V (t) d t

и в результате вычисляется одна точка ВАХ для заданного значения $I$ . Отметим, что при усреднении интегрирования начинается от времени $T m i n$ , после которого решение стабилизируется. Далее меняется значение тока на δ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 149 ms, sys: 633 µs, total: 150 ms
Wall time: 160 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()

_images/926157b9b1cbbf4e54f877fbf31dcb8674cfc2a62a2b5d418cd342dd3229fb31.png

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()

_images/6bdb9d2930f5ef4852208225d8159e3867cefe5daa55e9b516607dc4c0c495c6.png

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)

[ 5.85144937e-12 -8.31711194e-01  2.25505248e-13  8.31711194e-01
  2.06995079e-13]

print(res[0])

5.851449372528068e-12

print(res[1:])

[-8.31711194e-01  2.25505248e-13  8.31711194e-01  2.06995079e-13]

Вычислим ВАХ

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)

26.542141675949097

График ВАХ

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()

_images/4218258008c5e47029e0b3841099cb6d9f89bf547707ae470c1e77a5488b9690.png

Описание ВАХ#

Как видно из ВАХ с увеличением тока начиная от нуля до 0.8 значения напряжение равна нулю (стационарный эффект Джозефсона). При достижении значения тока 0.8 появляется напряжение и далее с ростом тока растет и напряжение. Как было отмечено в описании модели ток при котором появляется напряжение назывется критическим током. Теперь рассмотрим как критический ток меняется с изменением внешнего потока.

Исследование влияния $φ_{e x t}$ на ВАХ#

Проведем вычисления при различных значений $φ_{e x t}$ , например при $φ_{e x t} = 0$ , $φ_{e x t} = 0.25$ и $φ_{e x t} = 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)

5.828256845474243

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)

13.274491786956787

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)

26.894043922424316

График ВАХ#

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()

_images/d750aab7cb8de9213e467e0bc1c6d1fe66a6b9f23689eb634036bec11a11e43c.png

Как видно из графика ВАХ при $φ_{e x t} = 0$ критический ток равен $I_{c} = 2$ , при $φ_{e x t} = 0.25$ - $I_{c} = 1.58$ , а при $φ_{e x t} = 0.5$ - $I_{c} = 0.8$ , т.е. с изменением $φ_{e x t}$ меняется величина критического тока.

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)

Execution time 1.759821891784668 s
2.0009999999998906

Исследование зависимости критического тока от внешнего магнитного поля#

Для этого вводим следующие параметры:

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)

0.04

Создаем массив значений внешнего поля для которых нужно вычислить критическй ток

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

CPU times: user 53.7 s, sys: 71.3 ms, total: 53.8 s
Wall time: 53.8 s

print(Ic_array)

[2.009 1.993 1.949 1.883 1.8   1.704 1.598 1.486 1.368 1.248 1.124 1.
 0.877 0.877 1.    1.124 1.248 1.368 1.486 1.598 1.704 1.8   1.883 1.949
 1.993 2.009 1.993 1.949 1.883 1.8   1.704 1.598 1.486 1.368 1.248 1.124
 1.    0.877 0.877 1.    1.124 1.248 1.368 1.486 1.598 1.704 1.8   1.883
 1.949 1.993 2.009]

Посторим график зависимости критического тока от внешнего поля#

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()

_images/c3ca86122f5346d51da5812a6572fdd29e7fa943db2cd8570d63e4a3a877167d.png

Данный график демонстрирует модулирование критического тока с изменением $φ_{e x t}$ .

Параллелизация с 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

set_num_threads(5)

config.NUMBA_NUM_THREADS = 20

config.NUMBA_NUM_THREADS

config.THREADING_LAYER = 'omp'

threading_layer()

'tbb'

config.THREADING_LAYER

'omp'

# demonstrate the threading layer chosen
print("Threading layer chosen: %s" % threading_layer())

Threading layer chosen: tbb

get_num_threads()

set_num_threads(5)

get_num_threads()

import os
print('NUMBA NUM THREADS')
try:
    print(os.environ['NUMBA_NUM_THREADS'])
except:
    print(config.NUMBA_DEFAULT_NUM_THREADS)

NUMBA NUM THREADS
80

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)

CPU times: user 56.3 s, sys: 45 ms, total: 56.4 s
Wall time: 13.8 s

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()

_images/60c32fbb30e4846a05f4040a3fccd52e06c8f9768b960127caa0a5ad6fd6f757.png

Вычисления при разных чисел потоков.

На этом этапе реализована параллелизация вычислении критического тока. На графике представлена рассчитанная зависимость критического тока СКВИДа от потока $φ_{e x t}$ в интервале потока [0, 2]. Данная зависимость демонстрирует модулирование критического тока с изменением $φ_{e x t}$ , что согласуется с известными результатами. Теперь можно вычислить зависимость времени вычисления от числа вычислительных потоков.

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 52.91446542739868 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.294822692871094 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 8.151139974594116 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()

_images/2137997a657d947b07394fc2febf0a2456242d3ead1fce59204519e7f5a89365.png

Разработанный инструментарий позволяет моделировать физические свойства сверхпроводникового квантового интерферометра с двумя джозефсоновскими переходами (СКВИД постоянного тока). Реализованы алгоритмы для вычисления вольт-амперной характеристики СКВИДа под воздействием внешнего магнитного поля и зависимости критического тока СКВИДа от внешнего магнитного поля. При программной реализации разработанных алгоритмов использовались функции библиотеки Numba, в том числе механизм распараллеливания при вычислении зависимости критического тока от потока внешнего магнитного поля. Представленные результаты были получены на экосистеме ML/DL/HPC Гетерогенной платформы HybriLIT (ЛИТ ОИЯИ).

Работа выполнена при финансовой поддержки РНФ в рамках проекта No 22-71-10022.

Contents