.. _akustika_ravan:

Rešavanje na domenu oblika kvadrata
==========================================
 
Krenućemo od najjednostavnijeg dvodimenzionog slučaja stojećeg talasa u akustici. Za talasni broj :math:`k_0=2 \pi n` za :math:`n=2`, treba rešiti Helmholcovu (*Helmholtz*) jednačinu oblika:

.. math:: 
    -\frac{\partial^2 u}{\partial x^2} - \frac{\partial^2 u}{\partial y^2} - k_0^2 u = f, \quad \Omega=[0,1]^2,

uz Dirihleove granične uslove 

.. math:: 
    u(x,y)=0, \quad (x,y) \in \partial \Omega

prikazane na :numref:`akustika-2d` i član koji specificira izvor :math:`f(x,y)=k_0^2 \sin(k_0x) \sin(k_0 y)`.

.. _akustika-2d:

.. figure:: akustika1.png
    :width: 60%

    Postavka problema i granični uslovi.

Postoji analitičko rešenje ovog problema i ono glasi:

.. math:: 
    u(x,y) = \sin(k_0 x) \sin(k_0 y).

Za više detalja u pogledu teorijske pozadine diferencijalne jednačine i graničnih uslova čitalac može konsultovati :cite:t:`ihlenburg1998finite`. Rešenje metodom konačnih elemenata (MKE) je takođe dostupno u okviru `Dolfinx tutorijala <https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_helmholtz.py>`_. 


Implementacija
----------------

Rešenje direktnog problema prikazano je na sledećem listingu. 

.. code-block:: python
    :caption: Rešenje problema prostiranja stojećeg talasa u 2D korišćenjem DeepXDE biblioteke
    :linenos:

    import deepxde as dde
    import numpy as np

    # Frekvencija
    n = 2

    precision_train = 10
    precision_test = 30
    weights = 100
    iterations = 10000
    learning_rate, num_dense_layers, num_dense_nodes, activation = 1e-3, 3, 150, "sin"

    # Uvezi sinus
    from deepxde.backend import tf
    sin = tf.sin

    # Osnovna PDE
    def pde(x, u):
        du_xx = dde.grad.hessian(u, x, i=0, j=0)
        du_yy = dde.grad.hessian(u, x, i=1, j=1)

        f = k0 ** 2 * sin(k0 * x[:, 0:1]) * sin(k0 * x[:, 1:2])
        return -du_xx - du_yy - k0 ** 2 * u - f

    # Egzaktno resenje
    def func(x):
        return np.sin(k0 * x[:, 0:1]) * np.sin(k0 * x[:, 1:2])

    # Da li je kol. tacka na granici?
    def boundary(_, on_boundary):
        return on_boundary

    # Geometrija jedinicnog kvadrata
    geom = dde.geometry.Rectangle([0, 0], [1, 1])
    # Talasni broj
    k0 = 2 * np.pi * n
    # Talasna duzina
    wave_len = 1 / n

    hx_train = wave_len / precision_train
    nx_train = int(1 / hx_train)

    hx_test = wave_len / precision_test
    nx_test = int(1 / hx_test)

    # Dirihleov granicni uslov y=0 na granicama
    bc = dde.icbc.DirichletBC(geom, lambda x: 0, boundary)

    data = dde.data.PDE(
        geom,
        pde,
        bc,
        num_domain=nx_train ** 2,
        num_boundary=4 * nx_train,
        solution=func,
        num_test=nx_test ** 2,
    )

    # Mreza i model
    net = dde.nn.FNN([2] + [num_dense_nodes] * num_dense_layers + [1], activation, "Glorot uniform")
    model = dde.Model(data, net)

    # Forsiraj vece tezine za granicne uslove nego za unutrasnjost domena
    loss_weights = [1, weights]

    model.compile("adam", lr=learning_rate, metrics=["l2 relative error"], loss_weights=loss_weights)

    losshistory, train_state = model.train(iterations=iterations)
    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

Nakon standardnog importa odgovarajućih modula, počinjemo specifikacijom opštih parametara. Ovaj primer ima par specifičnosti u odnosu na ostale. Naime, da bi se uspešno modelovale talasne pojave pomoću NMPFZ, gustina kolokacionih tačaka mora da bude direktno proporcionalna frekvenciji. Što je viša frekvencija ``n``, manja je talasna dužina ``wave_len``, pa je potrebno više kolokacionih tačaka da pokrije domen. Ovde smo uzeli 10 kolokacionih tačaka po talasnoj dužini tokom treninga i 30 tačaka po talasnoj dužini u test skupu. 

.. code-block:: python

    # Frekvencija talasa
    n = 2
    precision_train = 10
    precision_test = 30
    weights = 100
    learning_rate, num_dense_layers, num_dense_nodes, activation = 1e-3, 3, 150, "sin"

Takođe, vidimo da koristimo arhitekturu sa manjim brojem slojeva, ali sa više neurona po sloju, kao i aktivacionu funkciju :math:`\sin(x)` koja bi trebalo da bude pogodnija za oponašanje talasnih fenomena. 

Sledi specifikacija same parcijalne diferencijalne jednačine u obliku funkcije gubitka kako smo to već navikli:

.. code-block:: python

    def pde(x, u):
        du_xx = dde.grad.hessian(u, x, i=0, j=0)
        du_yy = dde.grad.hessian(u, x, i=1, j=1)

        f = k0 ** 2 * sin(k0 * x[:, 0:1]) * sin(k0 * x[:, 1:2])
        return -du_xx - du_yy - k0 ** 2 * u - f

Ovde koristimo uslužnu funkciju ``dde.grad.hessian`` odabirom koordinate koja se diferencira i kojom se diferencira. U ovom primeru su granični uslovi elementarni, pa ih ovde nećemo posebno navoditi. 

Geometrija, talasni broj :math:`k_0=2 \pi \nu` i talasna dužina :math:`\lambda=\frac{1}{\nu}` daju se kao:

.. code-block:: python

    geom = dde.geometry.Rectangle([0, 0], [1, 1])
    k0 = 2 * np.pi * n
    wave_len = 1 / n

Jedina specifičnost koju dodatno treba naglasiti je da ponekad treba forsirati poštovanje graničnih uslova time što ćemo članu funkcije gubitka koji se odnosi na Dirihleov granični uslov dobiti veću težinu u odnosu na član koji se odnosi na diferencijalnu jednačinu. 

.. code-block:: python

    weights = 100
    loss_weights = [1, weights]
    model.compile("adam", lr=learning_rate, metrics=["l2 relative error"], loss_weights=loss_weights)

Da bi se izbegao ovaj korak koji sa sobom nosi eksperimentisanje sa različitim vrednostima težinskog faktora, granični uslov se kod *DeepXDE* može zadavati i direktnom transformacijom funkcije gubitka, ali ovde se time nećemo baviti.


Rezultati
----------------

Nakon 10.000 epoha obučavanja optimizacionom metodom ``Adam`` koji je protekao kao što je prikazano na :numref:`rezultati1-loss`, dobijamo stojeći talas čiji 3D prikaz možemo videti na :numref:`rezultati1`.

.. _rezultati1-loss:

.. figure:: rezultati1-loss.png
    :width: 80%

    Tok obučavanja NMPFZ

.. _rezultati1:

.. figure:: rezultati1.png
    :width: 80%

    Trodimenzioni prikaz talasa u domenu oblika kvadrata

Mera greške modela RMSE (*Root Mean Squared Error*) iznosi :math:`7,98 \cdot 10^{-2}`. Uz obraćanje posebne pažnje na forsiranje graničnih uslova, zatim arhitekturu NMPFZ i najzad tip aktivacione funkcije, uspeli smo da dobijemo prilično dobro rešenje. Čitalac može samostalno da proba kako bi promena frekvencije (a samim tim i talasne dužine), gustine kolokacionih tačaka, arhitekture, uticala na proces obučavanja modela. 

Ovde možemo dati i kratku preporuku **kako pristupiti modelovanju složenijih pojava**, sa složenijom geometrijom i kompleksnijim graničnim uslovima. Pošto NMPFZ rešavanje zavisi od većeg broja hiper-parametara, preporuka je da se prvo reši do kraja pojednostavljen problem baziran na istoj diferencijalnoj jednačini, ali sa jednostavnijom geometrijom i graničnim uslovima. Kada se stekne slika o tome koja kombinacija hiper-parametara vodi do konvergencije rešenja, onda je lakše pristupiti glavnom, kompleksnom problemu. S druge strane, postoji nekoliko alata koji pretragu hiper-parametara čine efikasnijom, kao već pomenuti `BlackFox <https://blackfox.ai>`_ koji koristi distribuirani genetski algoritam.