.. _akustika_socivo: Dvodimenzioni problem sa sočivom ====================================== Ovaj odeljak se bavi nešto modifikovanim `primerom iz tutorijala `_ za poznati softverski okvir za rad sa metodom konačnih elemenata **Deal.II**. Prvobitna svrha ovog primera je da simulira svojstva fokusiranja ultrazvučnog talasa koji generiše sočivo pretvarača sa promenljivom geometrijom. Savremene primene u medicinskom imidžingu koriste ultrazvučne talase ne samo za svrhe snimanja, već i za izazivanje određenih lokalnih efekata u materijalu, kao što su promene u optičkim svojstvima, koje se zatim mogu meriti drugim tehnikama snimanja. Vitalni sastojak ovih metoda je sposobnost fokusiranja intenziteta ultrazvučnog talasa u određenom delu materijala, idealno u tački, kako bi se mogla ispitati svojstva materijala na toj lokaciji. Karakteristična talasna dužina ultrazvuka nešto manja od fenomena koje smo do sada modelovali metodom NMPFZ, te bi nam bio potreban izuzetno veliki broj kolokacionih tačaka što produžava treniranje. Zato smo talasnu dužinu nešto povećali (smanjili frekvenciju), dok ostale parametre nismo menjali. Kako bismo izveli jednačine za ovaj problem, zvuk uzimamo kao talas kojim se širi promena pritiska: .. math:: \frac{\partial^2 U}{\partial t^2} - c^2 \Delta U = 0, gde je :math:`c` brzina zvuka, koja se zbog jednostavnosti uzima kao konstanta, :math:`U = U(x,t) ,\; x \in \Omega, \; t \in \mathrm{R}`. Postavka problema data je na :numref:`akustika-socivo`. .. _akustika-socivo: .. figure:: akustika3.png :width: 60% Postavka problema i granični uslovi Granica :math:`\Gamma=\partial\Omega` podeljena je na dva dela, i to :math:`\Gamma_1` i :math:`\Gamma_2 = \Gamma\setminus\Gamma_1`, gde :math:`\Gamma_1` predstavlja sočivo, a :math:`\Gamma_2` apsorbujuću granicu. Zapravo, želimo da napravimo takav granični uslov na :math:`\Gamma_2` tako da se oponaša znatno veći domen. Na :math:`\Gamma_1`, pretvarač generiše talase konstantne frekvencije :math:`\omega > 0` i konstantne jedinične amplitude: .. math:: U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1 Pošto nema drugih (internih ili graničnih) izvora i pošto samo izvor emituje talase frekvencije :math:`\omega`, dozvoljeno je da izvršimo razdvajanje promenljivih :math:`U(x,t) = \textrm{Re}\left(u(x)\,e^{i\omega t})\right)`. Kompleksna funkcija :math:`u(x)` opisuje prostornu zavisnost amplitude i faze (relativno u odnosu na izvor) talasa frekvencije :math:`\omega`, dok je amplituda veličina koja nas interesuje. Ako ovako formulisanu funkciju uvrstimo u talasnu jednačinu, vidimo da za :math:`u` imamo .. math:: -\omega^2 u(x) - c^2\Delta u(x) = 0, \qquad x \in \Omega, \\ u(x) = 1, \qquad x \in \Gamma_1. Da bismo našli odgovarajuće granične uslove na :math:`\Gamma_2` koji oponašaju apsorbujuću granicu, razmotrimo talas oblika :math:`V(x,t)=e^{i(k\cdot x - \omega t)}` frekvencije :math:`\omega` koji se prostire u pravcu :math:`k \in {\mathrm{R}^2}`. Da bi :math:`V` bio rešenje talasne jednačine, mora da važi :math:`|k|=\frac{\omega}{c}`. Pretpostavimo da talas dolazi do :math:`x_0 \in \Gamma_2` pod pravim uglom, na primer :math:`n=\frac{k}{|k|}` gde :math:`n` označava normalu na :math:`\Omega \in x_0`. Onda u :math:`x_0`, ovaj talas zadovoljava jednačinu .. math:: c (n\cdot\nabla V) + \frac{\partial V}{\partial t} = (i\, c\, |k| - i\, \omega) V = 0. Postavljanjem graničnog uslova .. math:: c (n\cdot\nabla U) + \frac{\partial U}{\partial t} = 0, \qquad x\in\Gamma_2, talasi koji udaraju u granicu :math:`\Gamma_2` pod pravim uglom biće savršeno apsorbovani. Sa druge strane, oni delovi talasnog polja koji ne padaju pod pravim uglom ne zadovoljavaju ovaj uslov, pa će dolaziti do parcijalnih refleksija. U osnovi, direktni delovi talasa će proći kroz granicu kao da ona ne postoji, dok će ostali biti reflektovani nazad u domen. Ukoliko smo spremni da prihvatimo ovako predloženu aproksimaciju, onda za :math:`u` važi sledeće: .. math:: :label: eq:sistem -\omega^2 u - c^2\Delta u = 0, \qquad x \in \Omega, \\ c (n\cdot\nabla u) + i\,\omega\, u=0, \qquad x \in \Gamma_2, \\ u=1, \qquad x \in \Gamma_1. prepoznajemo Helmholcovu jednačinu sa Dirihleovim uslovom na :math:`\Gamma_1` i mešanim graničnim uslovom na :math:`\Gamma_2`. Zbog uslova na :math:`\Gamma_2` ne možemo da tretiramo realne i imaginarne delove :math:`u` posebno. Ono što možemo da uradimo je da formiramo sistem od dve parcijalne diferencijalne jednačine u kojima figurišu realni i imaginarni deo :math:`u`, sa graničnim uslovima na :math:`\Gamma_2` koje vezuju ove dve komponente. Ako označimo da je :math:`v=\textrm{Re}\;u,\; w=\textrm{Im}\;u`, sistem :math:numref:`eq:sistem` glasi: .. math:: :label: eq:sistem1 -\omega^2 v - c^2 \Delta v=0 \quad \\ -\omega^2 w - c^2 \Delta w=0 \quad \\ x \in \Omega, \\ c (n\cdot\nabla v) - \omega \, w=0 \quad \\ c (n\cdot\nabla w) + \omega \, v=0 \quad \\ x \in \Gamma_2, \\ v=1 \quad \\ w=0 \quad \\ x \in \Gamma_1. Dakle, prve dve jednačine važe u celom domenu :math:`\Omega`, druge dve na granici :math:`\Gamma_2`, a poslednje dve na :math:`\Gamma_1`. Ovde prvi put imamo sistem diferencijalnih jednačina, ali ni to ne bi trebalo da bude problem za NMPFZ pristup, ako podrazumevamo da je sistem zatvoren, tj. jednoznačan. Implementacija ----------------- Na osnovu sistema jednačina :math:numref:`eq:sistem1` treba da formiramo kompozitnu funkciju gubitka, da formiramo NMPFZ mrežu i da je istreniramo na dovoljnom broju kolokacionih tačaka. Evo ključnih delova implementacije ostvarene pomoću okvira SCIANN: .. code-block:: python :caption: Rešenje problema prostiranja talasa u 2D domenu sa sočivom :linenos: import numpy as np import matplotlib.pyplot as plt import sciann as sn from numpy import pi from sciann.utils.math import diff, sign, sin, sqrt # Brzina talasa c = 1 # Frekvencija omega = 2*pi*4 x = sn.Variable('x') y = sn.Variable('y') v, w = sn.Functional (["v", "w"], [x, y], 3*[200] , 'sin') # Diferencijalne jednacine za v i w L1 = -omega**2 * v - c**2 * diff(v, x, order=2) - c**2 * diff(v, y, order=2) L2 = -omega**2 * w - c**2 * diff(w, x, order=2) - c**2 * diff(w, y, order=2) TOL = 0.015 # Dirihleov uslov na G1 (y=0 i 0.40.6) C9 = (1-sign(y - (0+TOL))) * ( (1 - sign(x-0.4)) + (1 + sign(x-0.6)) ) * ( -c*diff(v,y) - omega*w ) C10 = (1-sign(y - (0+TOL))) * ( (1 - sign(x-0.4)) + (1 + sign(x-0.6)) ) * ( -c*diff(w,y) + omega*v ) x_data, y_data = [], [] kolokacione_tacke = np.genfromtxt('kolokacione_tacke.txt', delimiter=" ") for e in kolokacione_tacke: ind, x1, y1 = e x_data.append(x1) y_data.append(y1) x_data, y_data = np.array(x_data), np.array(y_data) # Model i obucavanje m = sn.SciModel([x, y], [L1,L2,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10], 'mse', 'Adam') h = m.train([x_data, y_data], 12*['zero'], learning_rate=0.001, batch_size=1024, epochs=8000, adaptive_weights={'method':'NTK', 'freq':200}) # Test x_test, y_test = np.meshgrid( np.linspace(0, 1, 200), np.linspace(0, 1, 200) ) v_pred = v.eval(m, [x_test, y_test]) w_pred = w.eval(m, [x_test, y_test]) Nakon uobičajenih importa paketa, formiramo NMPFZ mreže za realni deo :math:`v` i imaginarni deo :math:`w` nepoznate funkcije :math:`u`. I ovde ćemo ići sa aktivacionom funkcijom :math:`\sin(x)`. Interesantan deo koda je definisanje graničnog uslova na :math:`\Gamma_1` prema poslednje dve jednačine u sistemu :math:numref:`eq:sistem1`: .. code-block:: python a,b,c = 0.39762422, -1.57715550, 1.60337246 C1 = (1 - sign(y - (a + b*x + c*x**2 + TOL))) * (1 + sign(x-0.4)) * (1 - sign(x-0.6)) * (1-v) C2 = (1 - sign(y - (a + b*x + c*x**2 + TOL))) * (1 + sign(x-0.4)) * (1 - sign(x-0.6)) * (w-0) Jednačina :math:`y = a + b \cdot x + c \cdot x^2` predstavlja jednačinu spoljne linije sočiva :math:`\Gamma_1`, u kojoj su koeficijenti :math:`a,b,c` dobijeni fitovanjem. Dakle, prva zagrada u graničnim uslovima znači da uzimamo kolokacione tačke koje pripadaju tankom pojasu iznad linije :math:`\Gamma_1`, dok druga i treća zagrada imaju nenultu vrednost samo ako je :math:`0,4 < x < 0,6`. Granični uslovi ``C3, C4, C5, C6, C7, C8`` se odnose na mešanu formulaciju prema druge dve jednačine u sistemu :math:numref:`eq:sistem1` i važe na :math:`\Gamma_2`. Na primer: .. code-block:: python C5 = (1+sign(x - (1-TOL))) * ( c*diff(v,x) - omega*w ) C6 = (1+sign(x - (1-TOL))) * ( c*diff(w,x) + omega*v ) se odnosi na desnu granicu gde je ``x=1`` i uzima kolokacione tačke koje se nalaze u tankom pojasu širine ``TOL`` sa leve strane te granice. Komponente funkcije gubitka ``C9`` i ``C10`` odnose se takođe na granicu :math:`\Gamma_2`, ali na liniji gde je :math:`y=0`` i :math:`x < 0,4` ili :math:`x > 0,6`: .. code-block:: python C9 = (1-sign(y - (0+TOL))) * ( (1 - sign(x-0.4)) + (1 + sign(x-0.6)) ) * ( -c*diff(v,y) - omega*w ) C10 = (1-sign(y - (0+TOL))) * ( (1 - sign(x-0.4)) + (1 + sign(x-0.6)) ) * ( -c*diff(w,y) + omega*v ) Ovim smo kompletirali svih 12 komponenti funkcije gubitka. Pošto ih ima toliko, nije jednostavno izvršiti njihovo ponderisanje, odnosno dodelu težina svakoj komponenti. U ovakvim situacijama pomažu metode za adaptivno određivanje težina komponenti tokom obuke. U našem rešenju: .. code-block:: python h = m.train([x_data, y_data], 12*['zero'], learning_rate=0.001, batch_size=1024, epochs=8000, adaptive_weights={'method':'NTK', 'freq':200}) upotrebili smo inovativnu metodu *Neural Tangent Kernel* (NTK) prema :cite:t:`wang2022and`. Objašnjenje metode izlazi iz okvira ovog praktikuma, pa je nećemo detaljno razrađivati. Takođe, valja napomenuti da smo kolokacione tačke učitali iz posebnog fajla ``kolokacione_tacke.txt``, koji je dobijen tako što smo ispisali čvorove konačnih elemenata koji se dobijaju iz generatora mreže `paketa Deal.II `_. Rezultati ------------- Kao što je na početku odeljka već rečeno, primer je preuzet iz dokumentacije za paket koji se bavi analizom metodom konačnih elemenata `Deal.II `_, tako da možemo da uporedimo rešenje za :math:`v=\textrm{Re}\;u` koje smo dobili pomoću NMPFZ (:numref:`socivo-res-pinn`) i rešenje koje se dobija klasičnom metodom konačnih elemenata (:numref:`socivo-res-fem`). Rešenje dobijeno MKE metodom možemo smatrati referentnim, jer je korišćena veoma gusta mreža i pokazana je konvergencija. .. _socivo-res-pinn: .. figure:: socivo-v.png :width: 80% Rešenje za :math:`v=\textrm{Re}\;u` dobijeno pomoću NMPFZ .. _socivo-res-fem: .. figure:: rezultati-socivo-fem1.png :width: 80% Rešenje za :math:`v=\textrm{Re} \; u` dobijeno metodom konačnih elemenata Kvalitativno gledano, NMPFZ rešenje ima iste karakteristike kao MKE rešenje. Međutim, oko sočiva i prima granicama :math:`\Gamma_2` očigledan je pad kvaliteta rešenja. Uzrok možemo tražiti na nekoliko mesta: - Prvo, kod MKE je izvor na liniji sočiva :math:`\Gamma_1` moguće preciznije specificirati po samoj liniji, a ne u pojedinačnim tačkama kao kod NMPFZ. - Moguće je da 50.000 kolokacionih tačaka nije dovoljno za obučavanje. - Primećeno je da obučavanje optimizacionim algoritam ``Adam`` ne može da spusti vrednost gubitka ispod neke granice. Ovde verovatno treba eksperimentisati sa varijabilnom stopom učenja, ili dodati neki drugi optimizator kao u primeru :ref:`poplavni`. Dalje eksperimentisanje na ovom primeru ostavljamo čitaocu.