מידול נוזל זורם ללא צמיגות באמצעות QUICK-PDE

הערה

Qiskit Functions הן תכונה ניסיונית הזמינה רק למשתמשי IBM Quantum® Premium Plan, Flex Plan ו-On-Prem (via IBM Quantum Platform API) Plan. הן נמצאות בשלב גרסת תצוגה מקדימה וכפופות לשינויים.

אומדן שימוש: 50 דקות על מעבד Heron r2. (הערה: זהו אומדן בלבד. זמן הריצה שלך עשוי להשתנות.)

שימו לב שזמן הביצוע של פונקציה זו הוא בדרך כלל יותר מ-20 דקות, לכן ייתכן שתרצו לפצל את המדריך הזה לשני חלקים: החלק הראשון, בו תקראו אותו ותפעילו את המשימות, והשני מספר שעות מאוחר יותר (תוך מתן זמן מספק למשימות להסתיים), כדי לעבוד עם תוצאות המשימות.

רקע

מדריך זה מבקש ללמד ברמה מבואית כיצד להשתמש בפונקציית QUICK-PDE לפתרון בעיות מולטי-פיזיקס מורכבות על יחידות עיבוד קוונטיות 156Q Heron R2 באמצעות H-DES (Hybrid Differential Equation Solver) של ColibriTD. האלגוריתם הבסיסי מתואר במאמר H-DES. שימו לב שפותר זה יכול גם לפתור משוואות לא-לינאריות.

בעיות מולטי-פיזיקס - כולל דינמיקת נוזלים, דיפוזיית חום ועיוות חומרים, בין היתר - ניתנות לתיאור נרחב באמצעות משוואות דיפרנציאליות חלקיות (PDEs).

בעיות כאלה רלוונטיות ביותר לתעשיות שונות ומהוות ענף חשוב במתמטיקה שימושית. עם זאת, פתרון משוואות דיפרנציאליות חלקיות מרובות משתנים לא-לינאריות מצומדות עם כלים קלאסיים נותר מאתגר בשל הדרישה לכמות משאבים גדולה באופן אקספוננציאלי.

פונקציה זו מתאימה למשוואות עם מורכבות ומשתנים הולכים וגדלים, והיא הצעד הראשון לפתיחת אפשרויות שבעבר נחשבו לבלתי פתירות. כדי לתאר במלואם בעיה שמידול שלה באמצעות PDEs, יש צורך לדעת את התנאים ההתחלתיים והגבוליים. אלה יכולים לשנות באופן ניכר את הפתרון של ה-PDE ואת הדרך למציאת הפתרון שלהם.

מדריך זה מלמד אותך כיצד:

להגדיר את הפרמטרים של פונקציית תנאי ההתחלה.
לכוונן את מספר הקיוביטים (המשמשים לקידוד הפונקציה של המשוואה הדיפרנציאלית), העומק ומספר הזריקות.
להריץ 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")

שלב 1: הגדרת מאפייני הבעיה לפתרון

מדריך זה מכסה את חוויית המשתמש משתי נקודות מבט: הבעיה הפיזיקלית הנקבעת על ידי התנאים ההתחלתיים, והרכיב האלגוריתמי בפתרון דוגמה של דינמיקת נוזלים על מחשב קוונטי.

לדינמיקת נוזלים חישובית (CFD) יש מגוון רחב של יישומים, ולכן חשוב ללמוד ולפתור את ה-PDEs הבסיסיים. משפחה חשובה של PDEs הן משוואות Navier-Stokes, שהן מערכת של משוואות דיפרנציאליות חלקיות לא-לינאריות המתארות את תנועת הנוזלים. הן רלוונטיות ביותר לבעיות מדעיות ויישומים הנדסיים.

בתנאים מסוימים, משוואות Navier-Stokes מצטמצמות למשוואת Burgers, משוואת הסעה-דיפוזיה המתארת תופעות המתרחשות בדינמיקת נוזלים, דינמיקת גזים ואקוסטיקה לא-לינארית, בין היתר, באמצעות מידול מערכות דיסיפטיביות.

הגרסה החד-ממדית של המשוואה תלויה בשני משתנים: $t \in \mathbb{R}_{\geq 0}$ המדגימה את הממד הזמני, $x \in \mathbb{R}$ המייצגת את הממד המרחבי. הצורה הכללית של המשוואה נקראת משוואת Burgers הצמיגה ונכתבת:

$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial^2 u}{\partial^2 x},$

כאשר $u(x,t)$ הוא שדה מהירות הנוזל במיקום נתון $x$ ובזמן $t$ , ו- $\nu$ הוא הצמיגות של הנוזל. צמיגות היא תכונה חשובה של נוזל המודדת את ההתנגדות התלויה בקצב שלו לתנועה או לעיוות, ולכן היא ממלאת תפקיד מכריע בקביעת הדינמיקה של נוזל. כאשר הצמיגות של הנוזל היא אפס ( $\nu$ = 0), המשוואה הופכת למשוואת שימור שיכולה לפתח אי-רציפויות (גלי הלם), עקב היעדר ההתנגדות הפנימית שלה. במקרה זה, המשוואה נקראת משוואת Burgers הבלתי-צמיגה והיא מקרה מיוחד של משוואת גל לא-לינארית.

למען הדיוק, זרימות בלתי-צמיגות אינן מתרחשות בטבע, אך כאשר ממדלים זרימה אווירודינמית, עקב ההשפעה האינפיניטסימלית של ההובלה, שימוש בתיאור בלתי-צמיג של הבעיה יכול להיות שימושי. באופן מפתיע, יותר מ- 70% מתיאוריית האווירודינמיקה עוסקת בזרימות בלתי-צמיגות.

מדריך זה משתמש במשוואת Burgers הבלתי-צמיגה כדוגמה CFD לפתרון על יחידות עיבוד קוונטיות של IBM® באמצעות QUICK-PDE, הניתנת על ידי המשוואה:

$\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0$

תנאי ההתחלה לבעיה זו מוגדר לפונקציה לינארית: $u(t=0,x) = ax + b,\text{ with }a,b\in\mathbb{R}$ כאשר $a$ ו- $b$ הם קבועים שרירותיים המשפיעים על צורת הפתרון. אתם יכולים לכוונן את $a$ ו- $b$ ולראות כיצד הם משפיעים על תהליך הפתרון ועל הפתרון.

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  ])}}

שלב 2 (במידת הצורך): אופטימיזציה של הבעיה לביצוע על חומרה קוונטית

כברירת מחדל, הפותר משתמש בפרמטרים מבוססי-פיזיקה, שהם פרמטרי מעגל התחלתיים עבור מספר קיוביטים ועומק נתונים שממנו הפותר שלנו יתחיל.

הזריקות הן גם חלק מהפרמטרים עם ערך ברירת מחדל, מכיוון שכיוונן עדין הוא חשוב.

בהתאם לתצורה שאתם מנסים לפתור, ייתכן שיהיה צורך להתאים את פרמטרי האלגוריתם כדי להשיג פתרונות מספקים; לדוגמה, הוא יכול לדרוש יותר או פחות קיוביטים למשתנה $t$ ו- $x$ , בהתאם ל- $a$ ו- $b$ . הדברים הבאים מכווננים את מספר הקיוביטים לפונקציה למשתנה, את העומק לפונקציה ואת מספר הזריקות.

אתם יכולים גם לראות כיצד לציין את ה-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  ])}}

שלב 3: השוואת ביצועי האלגוריתם

אתם יכולים להשוות את תהליך ההתכנסות של הפתרון שלנו (HDES) של job_2 לביצועים של אלגוריתם ופותר רשתות נוירונים מבוססי-פיזיקה (PINN) (ראו את המאמר ואת מאגר GitHub המשויך).

בדוגמה של הפלט של job_2 (גישה מבוססת-קוונטים), רק 13 פרמטרים (12 פרמטרי מעגל פלוס פרמטר קנה מידה אחד) מאופטמזים עם הפותר הקלאסי. תהליך ההתכנסות הוא כדלקמן:

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

זה אומר שניתן להגיע להפסד מתחת ל-0.0015 לאחר 28 איטרציות, ועם אופטימיזציה של מעט פרמטרים קלאסיים בלבד.

כעת נוכל להשוות את אותו הדבר לפתרון PINN עם התצורה המוגדרת כברירת המחדל שהוצעה על ידי המאמר באמצעות אופטימיזר מבוסס-גרדיאנט. המקבילה של המעגל שלנו עם 13 פרמטרים שיש לאפטמז היא הרשת הנוירונית, הדורשת לפחות שמונה שכבות של 20 נוירונים, ולכן כוללת אופטימיזציה של 3021 פרמטרים. לאחר מכן, ההפסד היעד מושג בשלב 315, הפסד: 0.0014988397.

גרף המציג נתוני PINN בהשוואה לפונקציית HDES-Qiskit.

כעת, מכיוון שאנו רוצים לעשות השוואה הוגנת, עלינו להשתמש באותו אופטימיזר בשני המקרים. המספר הנמוך ביותר של איטרציות שמצאנו עבור 12 שכבות של 20 נוירונים = 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]
   20 5.398521572351456e-02 1.0e+00 9.98e-03  1e-02  1e-02 0:02.3
   40 5.444650724530220e-02 1.0e+00 9.97e-03  1e-02  1e-02 0:05.1
   60 4.447407275438309e-02 1.0e+00 9.95e-03  1e-02  1e-02 0:08.2
   80 2.068969979882240e-02 1.0e+00 9.94e-03  1e-02  1e-02 0:11.7
  120 1.028892211616039e-02 1.0e+00 9.91e-03  1e-02  1e-02 0:20.1
  140 5.140972323715687e-03 1.0e+00 9.90e-03  1e-02  1e-02 0:25.4
  180 3.811701666563749e-03 1.0e+00 9.87e-03  1e-02  1e-02 0:37.4
  200 3.189878538250923e-03 1.0e+00 9.85e-03  1e-02  1e-02 0:44.2
  240 2.547040116041899e-03 1.0e+00 9.83e-03  1e-02  1e-02 0:59.7
  280 2.166548743844032e-03 1.0e+00 9.80e-03  1e-02  1e-02 1:18.0
  300 1.783065614290535e-03 1.0e+00 9.79e-03  1e-02  1e-02 1:28.4
  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())

שלב 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()

פלט של תא הקוד הקודם

סקר מדריך

אנא הקדישו דקה לספק משוב על מדריך זה. התובנות שלכם יעזרו לנו לשפר את הצעות התוכן שלנו וחוויית המשתמש:

קישור לסקר

רקע​

דרישות​

הגדרה​

שלב 1: הגדרת מאפייני הבעיה לפתרון​

שלב 2 (במידת הצורך): אופטימיזציה של הבעיה לביצוע על חומרה קוונטית​

שלב 3: השוואת ביצועי האלגוריתם​

שלב 4: שימוש בתוצאה​

סקר מדריך​

רקע