주 콘텐츠로 건너뛰기

QUICK-PDE를 사용하여 비점성 유체 흐름 모델링

참고

Qiskit Functions는 IBM Quantum® Premium Plan, Flex Plan, On-Prem(IBM Quantum Platform API를 통한) Plan 사용자에게만 제공되는 실험적 기능입니다. 현재 미리 보기 릴리스 상태이며 변경될 수 있습니다.

사용량 예상: Heron r2 프로세서 기준 50분. (참고: 이는 예상치일 뿐이며, 실제 실행 시간은 다를 수 있습니다.)

이 함수의 실행 시간은 일반적으로 20분을 초과하므로, 이 튜토리얼을 두 섹션으로 나누는 것을 권장합니다. 첫 번째 섹션에서는 내용을 읽고 작업을 실행하고, 두 번째 섹션은 몇 시간 후(작업이 완료될 충분한 시간을 두고)에 진행하여 작업 결과를 처리합니다.

배경

이 튜토리얼은 ColibriTD의 H-DES(하이브리드 미분 방정식 솔버)를 사용하여 156Q Heron R2 QPU에서 복잡한 다중 물리학 문제를 해결하는 QUICK-PDE 함수의 사용법을 입문 수준에서 설명합니다. 기반 알고리즘은 H-DES 논문에 설명되어 있습니다. 이 솔버는 비선형 방정식도 풀 수 있습니다.

유체 역학, 열 확산, 재료 변형 등을 포함한 다중 물리학 문제는 편미분 방정식(PDE)으로 보편적으로 기술할 수 있습니다.

이러한 문제들은 다양한 산업에서 매우 중요하며, 응용 수학의 핵심 분야를 구성합니다. 그러나 고전적인 도구로 비선형 다변수 연립 PDE를 푸는 것은 지수적으로 방대한 양의 자원이 필요하기 때문에 여전히 어렵습니다.

이 함수는 복잡성과 변수가 증가하는 방정식에 적합하며, 한때 풀 수 없다고 여겨졌던 문제들을 해결할 가능성을 여는 첫 번째 단계입니다. PDE로 모델링된 문제를 완전히 기술하려면 초기 조건과 경계 조건을 알아야 합니다. 이 조건들은 PDE의 해와 해를 찾는 경로를 크게 바꿀 수 있습니다.

이 튜토리얼에서 배울 내용:

  1. 초기 조건 함수의 매개변수를 정의하는 방법
  2. Qubit 수(미분 방정식의 함수를 인코딩하는 데 사용), 깊이, 샷 수를 조정하는 방법
  3. QUICK-PDE를 실행하여 기저 미분 방정식을 푸는 방법

요구 사항

이 튜토리얼을 시작하기 전에 다음 항목이 설치되어 있는지 확인하세요:

  • Qiskit SDK v2.0 이상 (pip install qiskit)
  • Qiskit Functions Catalog (pip install qiskit-ibm-catalog)
  • Matplotlib (pip install matplotlib)
  • QUICK-PDE 함수 접근 권한. 접근 요청 양식을 작성하세요.

설정

API 키를 사용하여 인증하고 다음과 같이 함수를 선택하세요:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit-ibm-catalog
import numpy as np
import matplotlib.pyplot as plt
from qiskit_ibm_catalog import QiskitFunctionsCatalog

catalog = QiskitFunctionsCatalog(
channel="ibm_quantum_platform",
instance="INSTANCE_CRN",
token="YOUR_API_KEY", # Use the 44-character API_KEY you created and saved from the IBM Quantum Platform Home dashboard
)

quick = catalog.load("colibritd/quick-pde")

Step 1: 풀어야 할 문제의 속성 설정

이 튜토리얼은 두 가지 관점에서 사용자 경험을 다룹니다. 첫 번째는 초기 조건에 의해 결정되는 물리적 문제이고, 두 번째는 양자 컴퓨터에서 유체 역학 예제를 푸는 알고리즘적 구성 요소입니다.

전산 유체 역학(CFD)은 광범위한 응용 분야를 가지고 있으므로, 기저 PDE를 연구하고 해결하는 것이 중요합니다. PDE의 중요한 계열 중 하나는 나비에-스토크스(Navier-Stokes) 방정식으로, 유체의 운동을 기술하는 비선형 편미분 방정식 시스템입니다. 이 방정식은 과학적 문제와 공학적 응용에서 매우 중요합니다.

특정 조건에서 나비에-스토크스 방정식은 버거스(Burgers') 방정식으로 단순화됩니다. 버거스 방정식은 소산 시스템을 모델링함으로써 유체 역학, 기체 역학, 비선형 음향학 등에서 발생하는 현상을 기술하는 대류-확산 방정식입니다.

이 방정식의 1차원 버전은 두 변수에 의존합니다. tR0t \in \mathbb{R}_{\geq 0}는 시간 차원을 모델링하고, xRx \in \mathbb{R}은 공간 차원을 나타냅니다. 방정식의 일반 형태는 점성 버거스 방정식이라고 하며 다음과 같습니다:

ut+uux=ν2u2x,\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial^2 x},

여기서 u(x,t)u(x,t)는 주어진 위치 xx와 시간 tt에서의 유체 속도 필드이고, ν\nu는 유체의 점도입니다. 점도는 운동 또는 변형에 대한 속도 의존적 저항을 측정하는 유체의 중요한 특성이며, 따라서 유체 역학의 결정에 중요한 역할을 합니다. 유체의 점도가 0(ν\nu = 0)이면, 방정식은 내부 저항의 부재로 인해 불연속(충격파)이 발생할 수 있는 보존 방정식이 됩니다. 이 경우, 방정식은 비점성 버거스 방정식이라 하며, 비선형 파동 방정식의 특수한 경우입니다.

엄밀히 말하면, 비점성 흐름은 자연에서 발생하지 않습니다. 그러나 공기역학적 흐름을 모델링할 때, 수송의 미소한 효과로 인해 비점성 설명을 사용하는 것이 유용할 수 있습니다. 놀랍게도, 공기역학 이론의 70% 이상이 비점성 흐름을 다룹니다.

이 튜토리얼에서는 QUICK-PDE를 사용하여 IBM® QPU에서 CFD 예제로 비점성 버거스 방정식을 풀며, 방정식은 다음과 같습니다:

ut+uux=0\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0

이 문제의 초기 조건은 선형 함수로 설정됩니다: u(t=0,x)=ax+b, with a,bRu(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R} 여기서 aabb는 해의 형태에 영향을 주는 임의의 상수입니다. aabb를 조정하여 풀이 과정과 해에 미치는 영향을 확인할 수 있습니다.

job = quick.run(
use_case="cfd",
physical_parameters={"a": 1.0, "b": 1.0},
)
print(job.result())
{'functions': {'u': array([[1.        , 0.96112378, 0.9230742 , 0.88616096, 0.85058445,
0.81644741, 0.78376878, 0.75249908, 0.72253689, 0.69374562,
0.66597013, 0.63905258, 0.61284684, 0.58723093, 0.56211691,
0.53745752, 0.51324915, 0.48953036, 0.46637547, 0.44388257,
0.4221554 , 0.40127848, 0.38128488, 0.36211604, 0.34357308,
0.32525895, 0.30651089, 0.28632252, 0.26325504, 0.23533692],
[1.2375 , 1.19267729, 1.14850734, 1.10544526, 1.06382155,
1.02385326, 0.98565757, 0.94926734, 0.91464784, 0.88171402,
0.85034771, 0.82041411, 0.79177677, 0.76431068, 0.73791248,
0.71250742, 0.68805224, 0.66453346, 0.64196021, 0.62035121,
0.59971506, 0.5800232 , 0.56117499, 0.54295419, 0.52497612,
0.50662498, 0.48698059, 0.4647339 , 0.43809065, 0.40466247],
[1.475 , 1.4242308 , 1.37394048, 1.32472956, 1.27705866,
1.23125911, 1.18754636, 1.1460356 , 1.10675879, 1.06968242,
1.03472529, 1.00177563, 0.9707067 , 0.94139043, 0.91370806,
0.88755732, 0.86285533, 0.83953655, 0.81754494, 0.79681986,
0.77727473, 0.75876792, 0.74106511, 0.72379234, 0.70637915,
0.687991 , 0.66745028, 0.64314527, 0.61292625, 0.57398802],
[1.7125 , 1.65578431, 1.59937362, 1.54401386, 1.49029576,
1.43866495, 1.38943515, 1.34280386, 1.29886974, 1.25765082,
1.21910288, 1.18313715, 1.14963664, 1.11847019, 1.08950364,
1.06260722, 1.03765842, 1.01453964, 0.99312968, 0.97328851,
0.95483439, 0.93751264, 0.92095522, 0.90463049, 0.88778219,
0.86935702, 0.84791997, 0.82155665, 0.78776186, 0.74331358],
[1.95 , 1.88733782, 1.82480676, 1.76329816, 1.70353287,
1.6460708 , 1.59132394, 1.53957212, 1.49098069, 1.44561922,
1.40348046, 1.36449867, 1.32856657, 1.29554994, 1.26529921,
1.23765712, 1.21246152, 1.18954273, 1.16871442, 1.14975716,
1.13239406, 1.11625736, 1.10084533, 1.08546864, 1.06918523,
1.05072304, 1.02838966, 0.99996803, 0.96259746, 0.91263913]])}, 'samples': {'t': array([0. , 0.03275862, 0.06551724, 0.09827586, 0.13103448,
0.1637931 , 0.19655172, 0.22931034, 0.26206897, 0.29482759,
0.32758621, 0.36034483, 0.39310345, 0.42586207, 0.45862069,
0.49137931, 0.52413793, 0.55689655, 0.58965517, 0.62241379,
0.65517241, 0.68793103, 0.72068966, 0.75344828, 0.7862069 ,
0.81896552, 0.85172414, 0.88448276, 0.91724138, 0.95 ]), 'x': array([0. , 0.2375, 0.475 , 0.7125, 0.95 ])}}

Step 2 (필요한 경우): 양자 하드웨어 실행을 위한 문제 최적화

기본적으로, 솔버는 물리적으로 정보화된 매개변수를 사용합니다. 이는 주어진 Qubit 수와 깊이에 대한 초기 Circuit 매개변수로, 솔버가 여기서부터 시작합니다.

샷도 기본값이 있는 매개변수의 일부이며, 세밀한 조정이 중요합니다.

풀려는 구성에 따라, 만족스러운 해를 얻기 위한 알고리즘의 매개변수를 조정해야 할 수 있습니다. 예를 들어, aabb에 따라 변수 ttxx에 대한 함수당 더 많거나 더 적은 Qubit이 필요할 수 있습니다. 다음은 함수당 변수별 Qubit 수, 함수당 깊이, 샷 수를 조정하는 방법을 보여줍니다.

Backend와 실행 모드를 지정하는 방법도 확인할 수 있습니다.

또한, 물리적으로 정보화된 매개변수가 최적화 과정을 잘못된 방향으로 이끌 수 있습니다. 이 경우, initialization 전략을 "RANDOM"으로 설정하여 비활성화할 수 있습니다.

job_2 = quick.run(
use_case="cfd",
physical_parameters={"a": 0.5, "b": 0.25},
nb_qubits={"u": {"t": 2, "x": 1}},
depth={"u": 3},
shots=[500, 2500, 5000, 10000],
initialization="RANDOM",
backend="ibm_kingston",
mode="session",
)
print(job_2.result())
{'functions': {'u': array([[0.25      , 0.24856543, 0.24687708, 0.2449444 , 0.24277686,
0.24038389, 0.23777496, 0.23495952, 0.23194702, 0.22874691,
0.22536866, 0.22182171, 0.21811551, 0.21425952, 0.2102632 ,
0.20613599, 0.20188736, 0.19752675, 0.19306361, 0.18850741,
0.18386759, 0.1791536 , 0.17437491, 0.16954096, 0.16466122,
0.15974512, 0.15480213, 0.1498417 , 0.14487328, 0.13990632],
[0.36875 , 0.36681313, 0.36457201, 0.36203594, 0.35921422,
0.35611615, 0.35275103, 0.34912817, 0.34525687, 0.34114643,
0.33680614, 0.33224532, 0.32747327, 0.32249928, 0.31733266,
0.31198271, 0.30645873, 0.30077002, 0.29492589, 0.28893564,
0.28280857, 0.27655397, 0.27018116, 0.26369944, 0.2571181 ,
0.25044645, 0.24369378, 0.23686941, 0.22998264, 0.22304275],
[0.4875 , 0.48506084, 0.48226695, 0.47912748, 0.47565158,
0.47184841, 0.46772711, 0.46329683, 0.45856672, 0.45354594,
0.44824363, 0.44266894, 0.43683103, 0.43073904, 0.42440212,
0.41782942, 0.4110301 , 0.4040133 , 0.39678818, 0.38936388,
0.38174955, 0.37395435, 0.36598742, 0.35785791, 0.34957498,
0.34114777, 0.33258544, 0.32389713, 0.315092 , 0.30617919],
[0.60625 , 0.60330854, 0.59996188, 0.59621902, 0.59208895,
0.58758067, 0.58270318, 0.57746549, 0.57187658, 0.56594545,
0.55968112, 0.55309256, 0.54618879, 0.53897879, 0.53147158,
0.52367614, 0.51560147, 0.50725658, 0.49865046, 0.48979211,
0.48069053, 0.47135472, 0.46179367, 0.45201638, 0.44203186,
0.4318491 , 0.42147709, 0.41092485, 0.40020136, 0.38931562],
[0.725 , 0.72155625, 0.71765682, 0.71331056, 0.70852631,
0.70331293, 0.69767926, 0.69163414, 0.68518643, 0.67834497,
0.6711186 , 0.66351618, 0.65554655, 0.64721855, 0.63854104,
0.62952285, 0.62017284, 0.61049986, 0.60051274, 0.59022035,
0.57963151, 0.56875509, 0.55759992, 0.54617486, 0.53448874,
0.52255042, 0.51036875, 0.49795257, 0.48531072, 0.47245205]])}, 'samples': {'t': array([0. , 0.03275862, 0.06551724, 0.09827586, 0.13103448,
0.1637931 , 0.19655172, 0.22931034, 0.26206897, 0.29482759,
0.32758621, 0.36034483, 0.39310345, 0.42586207, 0.45862069,
0.49137931, 0.52413793, 0.55689655, 0.58965517, 0.62241379,
0.65517241, 0.68793103, 0.72068966, 0.75344828, 0.7862069 ,
0.81896552, 0.85172414, 0.88448276, 0.91724138, 0.95 ]), 'x': array([0. , 0.2375, 0.475 , 0.7125, 0.95 ])}}

Step 3: 알고리즘 성능 비교

job_2의 솔루션(HDES)의 수렴 과정을 물리 정보 신경망(PINN) 알고리즘 및 솔버의 성능과 비교할 수 있습니다(논문 및 관련 GitHub 저장소 참조).

job_2의 출력(양자 기반 접근 방식) 예제에서는 클래식 솔버로 13개의 파라미터(Circuit 파라미터 12개 + 스케일링 파라미터 1개)만 최적화합니다. 수렴 과정은 다음과 같습니다:

optimizers:
CMA: {'ftarget': np.float64(0.1), 'verb_disp': 10, 'maxiter': 100}
CMA: {'ftarget': np.float64(0.005), 'verb_disp': 10, 'maxiter': 20}
CMA: {'ftarget': np.float64(0.0025), 'verb_disp': 10, 'maxiter': 30}
CMA: {'ftarget': np.float64(0.0005), 'verb_disp': 10, 'maxiter': 10}

500 shots
================== CMA =================
option: {'ftarget': np.float64(0.1), 'verb_disp': 10, 'maxiter': 100}
0/100, loss: 0.02456641

1000 shots
================== CMA =================
option: {'ftarget': np.float64(0.005), 'verb_disp': 10, 'maxiter': 20}
0/20, loss: 0.03641833
1/20, loss: 0.02461719
2/20, loss: 0.0283689
3/20, loss: 0.009898383
4/20, loss: 0.04454522
5/20, loss: 0.007019971
6/20, loss: 0.00811147
7/20, loss: 0.01592619
8/20, loss: 0.00764708
9/20, loss: 0.01401516
10/20, loss: 0.01767467
11/20, loss: 0.01220387

5000 shots
================== CMA =================
option: {'ftarget': np.float64(0.0025), 'verb_disp': 10, 'maxiter': 30}
0/30, loss: 0.01024792
1/30, loss: 0.004343748
2/30, loss: 0.01450951
3/30, loss: 0.008591284
4/30, loss: 0.00266414
5/30, loss: 0.007923613
6/30, loss: 0.02023853
7/30, loss: 0.01031438
8/30, loss: 0.009513116
9/30, loss: 0.008132266
10/30, loss: 0.005787766
11/30, loss: 0.00390582

10000 shots
================== CMA =================
option: {'ftarget': np.float64(0.0005), 'verb_disp': 10, 'maxiter': 10}
0/10, loss: 0.002386168
1/10, loss: 0.004024823
2/10, loss: 0.001311999
3/10, loss: 0.003433991
4/10, loss: 0.002339664
5/10, loss: 0.002978438
6/10, loss: 0.005458391
7/10, loss: 0.002026701
8/10, loss: 0.00207467
9/10, loss: 0.001947627
final_loss: 0.00151994463476429

즉, 28번의 반복만으로, 그리고 소수의 클래식 파라미터만 최적화하여 손실을 0.0015 이하로 낮출 수 있습니다.

이제 논문에서 제안된 기본 설정을 사용하는 PINN 솔루션(그래디언트 기반 옵티마이저 사용)과 동일한 결과를 비교할 수 있습니다. 13개의 파라미터를 최적화하는 우리의 Circuit에 해당하는 신경망은 최소 8개 레이어에 뉴런 20개가 필요하므로, 총 3021개의 파라미터를 최적화해야 합니다. 이 경우 목표 손실은 Step 315, loss: 0.0014988397에서 달성됩니다.

PINN 데이터와 HDES-Qiskit 함수를 비교한 그래프.

공정한 비교를 위해 두 경우 모두 동일한 옵티마이저를 사용해야 합니다. 뉴런 20개를 가진 12개 레이어 = 파라미터 4701개에서 찾아낸 최소 반복 횟수는 다음과 같습니다:

(10_w,20)-aCMA-ES (mu_w=5.9,w_1=27%) in dimension 4701 (seed=351961)
Iterat #Fevals function value axis ratio sigma min&max std t[m:s]
1 20 5.398521572351456e-02 1.0e+00 9.98e-03 1e-02 1e-02 0:02.3
2 40 5.444650724530220e-02 1.0e+00 9.97e-03 1e-02 1e-02 0:05.1
3 60 4.447407275438309e-02 1.0e+00 9.95e-03 1e-02 1e-02 0:08.2
4 80 2.068969979882240e-02 1.0e+00 9.94e-03 1e-02 1e-02 0:11.7
6 120 1.028892211616039e-02 1.0e+00 9.91e-03 1e-02 1e-02 0:20.1
7 140 5.140972323715687e-03 1.0e+00 9.90e-03 1e-02 1e-02 0:25.4
9 180 3.811701666563749e-03 1.0e+00 9.87e-03 1e-02 1e-02 0:37.4
10 200 3.189878538250923e-03 1.0e+00 9.85e-03 1e-02 1e-02 0:44.2
12 240 2.547040116041899e-03 1.0e+00 9.83e-03 1e-02 1e-02 0:59.7
14 280 2.166548743844032e-03 1.0e+00 9.80e-03 1e-02 1e-02 1:18.0
15 300 1.783065614290535e-03 1.0e+00 9.79e-03 1e-02 1e-02 1:28.4
16 320 2.045844215899706e-03 1.0e+00 9.78e-03 1e-02 1e-02 1:39.8
Stopping early: loss 0.001405 <= target 0.0015
CMA-ES finished. Best loss: 0.001404788694344461

job_2의 데이터로도 동일한 작업을 수행하여 PINN 솔루션과의 비교 그래프를 그릴 수 있습니다.

# check the loss function and compare between the two approaches
print(job_2.logs())

Step 4: 결과 활용

솔루션을 얻었으면 이제 원하는 방식으로 활용할 수 있습니다. 아래에서는 결과를 시각화하는 방법을 보여줍니다.

solution = job.result()

# Plot the solution of the second simulation job_2
_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution["samples"]["t"], solution["samples"]["x"])
ax.plot_surface(
t,
x,
solution["functions"]["u"],
edgecolor="royalblue",
lw=0.25,
rstride=26,
cstride=26,
alpha=0.3,
)
ax.scatter(t, x, solution, marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

이전 코드 셀의 출력

두 번째 실행에서의 초기 조건 차이와 그것이 결과에 미치는 영향을 확인해 보세요:

solution_2 = job_2.result()

# Plot the solution of the second simulation job_2
_ = plt.figure()
ax = plt.axes(projection="3d")

# plot the solution using the 3d plotting capabilities of pyplot
t, x = np.meshgrid(solution_2["samples"]["t"], solution_2["samples"]["x"])
ax.plot_surface(
t,
x,
solution_2["functions"]["u"],
edgecolor="royalblue",
lw=0.25,
rstride=26,
cstride=26,
alpha=0.3,
)
ax.scatter(t, x, solution_2, marker=".")
ax.set(xlabel="t", ylabel="x", zlabel="u(t,x)")

plt.show()

이전 코드 셀의 출력

튜토리얼 설문

이 튜토리얼에 대한 피드백을 남겨 주시면 감사하겠습니다. 여러분의 의견은 콘텐츠와 사용자 경험을 개선하는 데 큰 도움이 됩니다:

설문 링크

Note: This survey is provided by IBM Quantum and relates to the original English content. To give feedback on doQumentation's website, translations, or code execution, please open a GitHub issue.