Τι είναι η μέθοδος FDTD;

1. Εισαγωγή

Η μέθοδος πεπερασμένων διαφορών στο πεδίο του χρόνου (FDTD) επιλύει τις εξισώσεις Maxwell – οι οποίες περιγράφουν το φως – απευθείας στο πεδίο του χρόνου σε ένα διακριτό χωρικό πλέγμα. Σε αντίθεση με την ιχνηλάτηση ακτίνων, τους επιλυτές τρόπων ή τις προσεγγίσεις της μεθόδου της μεταφορικής μήτρας, η FDTD δεν κάνει απλοποιητικές υποθέσεις σχετικά με τη διάδοση ή τη δομή των τρόπων. Αντίθετα, τα ηλεκτρικά και μαγνητικά πεδία εξελίσσονται ρητά στο χρόνο μέσω της διακριτοποίησης των εξισώσεων στροβιλισμού – του νόμου του Ampère και του νόμου του Faraday – έτσι ώστε τα ηλεκτρομαγνητικά κύματα να διαδίδονται, να διασκορπίζονται, να παρεμβάλλονται και να διαθλούν με φυσικό τρόπο σε όλο το πεδίο προσομοίωσης. Αυτή η γενικότητα καθιστά την FDTD εξαιρετικά ακριβή και ευέλικτη, αλλά και υπολογιστικά απαιτητική όσον αφορά τη μνήμη και τον χρόνο εκτέλεσης.

\[ \nabla \times \boldsymbol{H} = \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t}, \qquad \nabla \times \boldsymbol{E} = - \mu \frac{\partial \boldsymbol{H}}{\partial t}. \]

Επειδή η μέθοδος FDTD δεν βασίζεται σε καμία υπόθεση σχετικά με την κατεύθυνση διάδοσης, τη δομή της συσκευής ή τη συμπεριφορά του πεδίου, μια μόνο προσομοίωση μπορεί να αποτυπώσει την απόκριση σε ευρύ φάσμα συχνοτήτων, τη διάθλαση και την συμβολή, τις επιπτώσεις του κοντινού πεδίου, τη γεωμετρία σε υπο-μήκη κύματος και τα παροδικά ηλεκτρομαγνητικά φαινόμενα, όλα εντός ενός ενιαίου πλαισίου. Αυτή η γενικότητα είναι αυτό που καθιστά την FDTD τόσο ισχυρή, αλλά έχει και ένα κόστος: το χωρικό πλέγμα πρέπει να επιλύει την μικρότερη φυσική κλίμακα μήκους, το χρονικό βήμα περιορίζεται αυστηρά από τη αριθμητική σταθερότητα και οι τρισδιάστατες προσομοιώσεις γίνονται γρήγορα απαιτητικές τόσο σε μνήμη όσο και σε χρόνο εκτέλεσης, ακόμη και για γεωμετρίες που φαίνονται μέτριες σε μέγεθος.

2. Γιατί δεν θα πρέπει να χρησιμοποιείτε την μέθοδο FDTD...

Στα περισσότερα πρακτικά οπτικά προβλήματα, η μέθοδος FDTD δεν θα πρέπει να είναι η πρώτη που θα χρησιμοποιήσετε. Απαιτεί σημαντικούς υπολογιστικούς πόρους και για πολλά συστήματα υπάρχουν πολύ πιο αποδοτικές προσεγγίσεις. Ο σχεδιασμός οπτικών συστημάτων, όπως κάμερες ή τηλεσκόπια, η μοντελοποίηση μεγάλης κλίμακας οπτικών διατάξεων ή η προσομοίωση της διαδρομής του φωτός σε μακροσκοπικές οπτικές διατάξεις αντιμετωπίζονται πολύ καλύτερα με τη χρήση της μεθόδου ιχνηλάτησης ακτίνων. Ομοίως, οι επίπεδες πολυεπίπεδες δομές, όπως τα λεπτού υμενίου ηλιακά κύτταρα, οι στοίβες OLED και οι οπτικές επιστρώσεις, αντιμετωπίζονται πολύ πιο αποδοτικά με τη χρήση μεθόδων μήτρας μεταφοράς. Σε αυτές τις περιπτώσεις, η FDTD προσθέτει σημαντικό υπολογιστικό κόστος χωρίς να παρέχει επιπλέον φυσική πληροφορία.

Οι περισσότερες μέθοδοι ηλεκτρομαγνητικής προσομοίωσης έχουν ως αφετηρία τις εξισώσεις του Maxwell. Τεχνικές όπως οι επιλυτές τρόπων και οι μέθοδοι μήτρας μεταφοράς προκύπτουν με την εισαγωγή απλοποιητικών υποθέσεων σχετικά με τη γεωμετρία, τη διάσταση ή τη δομή του πεδίου, με σκοπό την επιτάχυνση των υπολογισμών. Αυτές οι προσεγγίσεις είναι εξαιρετικά αποτελεσματικές, αλλά μόνο όταν οι υποκείμενες υποθέσεις τους ισχύουν.

Η μέθοδος FDTD βρίσκεται στην αντίθετη άκρη αυτού του φάσματος. Είναι εξαιρετικά ισχυρή, ιδιαίτερα για πολύπλοκα τρισδιάστατα συστήματα, ακανόνιστα γεωμετρικά σχήματα και προβλήματα όπου το πεδίο δεν μπορεί να αποσυντεθεί σε απλά modes ή στρώματα – όπως οι φωτονικές κρύσταλλοι, οι υποκυματικές δομές, τα μέσα που διασκορπίζουν έντονα και η σύζευξη κοντινού πεδίου. Το τίμημα για αυτή την ευελιξία είναι το υπολογιστικό κόστος, τόσο σε μνήμη όσο και σε χρόνο εκτέλεσης.

Το κλειδί είναι να επιλέξετε το κατάλληλο εργαλείο. Χρησιμοποιήστε απλοποιημένα μοντέλα όταν η φυσική το επιτρέπει και χρησιμοποιήστε τη μέθοδο FDTD όταν αυτές οι παραδοχές δεν ισχύουν πλέον. Το πλεονέκτημά της είναι η ευελιξία και η πληρότητα από φυσικής άποψης· το μειονέκτημα είναι το υπολογιστικό κόστος.

3. Οι εξισώσεις του Maxwell ως βάση για τη μέθοδο FDTD.

Στην ουσία του, το FDTD είναι πολύ απλό. Παίρνετε τις εξισώσεις του Maxwell, τις αναλύετε στις συνιστώσες τους για κάθε χωρική κατεύθυνση και στη συνέχεια επιλύετε αυτές τις εξισώσεις απευθείας σε ένα πλέγμα. Δεν υπάρχει κάποια μαγική διαδικασία – η μέθοδος είναι εννοιολογικά απλή και, σε πολλά σημεία, πιο απλή από πιο σύνθετες αριθμητικές προσεγγίσεις, όπως οι επιλυτές drift-diffusion. Για να ξεκινήσουμε, αναλύουμε τον νόμο του Ampère. Σε ένα γραμμικό μέσο με αγωγιμότητα \(\sigma\), διηλεκτρική σταθερά \(\epsilon\) και μαγνητική διαπερατότητα \(\mu\), αυτό μπορεί να γραφτεί ως εξής:

\[ \sigma \boldsymbol{E} + \epsilon \frac{\partial \boldsymbol{E}}{\partial t} = \nabla \times \boldsymbol{H} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ H_x & H_y & H_z \end{vmatrix} \]

το οποίο μπορεί να επεκταθεί ως

Ο νόμος του Φαραντέι διατυπώνεται ως εξής:

\[ -\sigma_m \boldsymbol{H} - \mu \frac{\partial \boldsymbol{H}}{\partial t} = \nabla \times \boldsymbol{E} = \begin{vmatrix} \hat{\boldsymbol{x}} & \hat{\boldsymbol{y}} & \hat{\boldsymbol{z}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ E_x & E_y & E_z \end{vmatrix} \]

το οποίο μπορεί να επεκταθεί ώστε να προκύψει:

Οι εξισώσεις που επισημαίνονται με μπλε και κόκκινο χρώμα αποτελούν τον ορισμό της μεθόδου FDTD και είναι η καρδιά της. Δεν υπάρχει κάτι πιο περίπλοκο από αυτό. Όλα τα υπόλοιπα είναι θέματα υλοποίησης: αντικατάσταση των παραγώγων με πεπερασμένες διαφορές, επιλογή οριακών συνθηκών, ορισμός πηγών φωτός (διέγερσης), εξαγωγή αποτελεσμάτων με τη χρήση ανιχνευτών και συγγραφή αποδοτικού (συχνά παράλληλου) κώδικα για την εξέλιξη των πεδίων στο χρόνο. Αυτά τα πρόσθετα στοιχεία καθιστούν έναν πρακτικό επιλυτή FDTD πιο περίπλοκο, αλλά δεν αλλάζουν την υποκείμενη απλότητα της βασικής μεθόδου.

4. Ένα απλό παράδειγμα FDTD για παιχνίδια.

Για να καταδείξουμε πόσο απλή είναι η μέθοδος FDTD – τόσο σε επίπεδο θεωρητικής κατανόησης όσο και στην πρακτική εφαρμογή – θα απλοποιήσουμε τις πλήρεις εξισώσεις σε ένα μονοδιάστατο, απλοποιημένο μοντέλο, το οποίο μπορεί να υλοποιηθεί σε λίγες δεκάδες γραμμές κώδικα σε MATLAB ή Python. Στόχος εδώ δεν είναι να αλλάξουμε τη φυσική, αλλά να αφαιρέσουμε την περιττή πολυπλοκότητα, ώστε η δομή της μεθόδου να είναι απόλυτα διαφανής.

Για να το επιτύχουμε αυτό, εξετάζουμε την μεταβολή μόνο κατά την κατεύθυνση \(x\) και ορίζουμε ρητά τις υπόλοιπες χωρικές παραγώγους σε μηδέν: \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. Αυτό μας αφήνει με ένα διαμήκες ηλεκτρομαγνητικό σύστημα με ηλεκτρικό πεδίο \(E_z(x,t)\) και μαγνητικό πεδίο \(H_y(x,t)\).

Υπό αυτές τις παραδοχές, οι εξισώσεις στροβιλισμού του Maxwell απλοποιούνται σε:

\[ \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} \]

\[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x}. \]

Αυτές είναι οι βασικές συνεχείς εξισώσεις για την μονοδιάστατη περίπτωση. Ωστόσο, σε αυτή τη μορφή δεν μπορούν να λυθούν άμεσα σε έναν υπολογιστή. Για να το κάνουμε αυτό, αντικαθιστούμε τις παραγώγους με πεπερασμένες διαφορές σε ένα διακριτό πλέγμα. Αυτό είναι η ουσία της μεθόδου: δείτε το πεπερασμένες διαφορές για μια γενική επισκόπηση.

Εισάγοντας ένα χωρικό πλέγμα \(x_i = i\Delta x\) και χρονικά βήματα \(t^n = n\Delta t\), και χρησιμοποιώντας την τυπική διαδοχική (leapfrog) διάταξη, οι εξισώσεις ενημέρωσης γίνονται:

\[ E_z^{n+1}(i) = E_z^n(i) + \frac{\Delta t}{\epsilon(i)\,\Delta x} \left[ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2}\right) - H_y^{n+\frac{1}{2}}\!\left(i-\frac{1}{2}\right) \right] \]

\[ H_y^{n+\frac{1}{2}}\!\left(i+\frac{1}{2}\right) = H_y^{n-\frac{1}{2}}\!\left(i+\frac{1}{2}\right) + \frac{\Delta t}{\mu\,\Delta x} \left[ E_z^n(i+1) - E_z^n(i) \right]. \]

Αυτό είναι ήδη η πλήρης μέθοδος FDTD στην απλούστερη μορφή της. Το ηλεκτρικό πεδίο ενημερώνεται με βάση τις χωρικές διαφορές στο μαγνητικό πεδίο και το μαγνητικό πεδίο ενημερώνεται με βάση τις χωρικές διαφορές στο ηλεκτρικό πεδίο, με τα δύο πεδία να εναλλάσσονται χρονικά.

Για ένα μη μαγνητικό μέσο, \(\mu = \mu_0\), η διηλεκτρική σταθερά σχετίζεται με τον δείκτη διάθλασης μέσω της εξίσωσης:

\[ \epsilon(x) = \epsilon_0 n(x)^2. \]

Αυτό σημαίνει ότι μια διαφορά στον δείκτη διάθλασης εισάγεται απλά αλλάζοντας την τιμή \(\epsilon\) σε όλη την περιοχή. Στο παρακάτω παράδειγμα, εξετάζουμε ένα μέσο που μεταβαίνει από την τιμή \(n=1\) στην τιμή \(n=3\), έτσι ώστε ένας παλμός που εκπέμπεται από την αριστερή πλευρά να ανακλάται εν μέρει και να διαπερνά εν μέρει στην επιφάνεια.

Στην πράξη, οι παραπάνω εξισώσεις απλώς τοποθετούνται μέσα σε έναν βρόχο χρονικής εξέλιξης: το μαγνητικό πεδίο ενημερώνεται με βάση το ηλεκτρικό πεδίο, εφαρμόζονται οι οριακές συνθήκες, στη συνέχεια το ηλεκτρικό πεδίο ενημερώνεται με βάση το μαγνητικό πεδίο και η διαδικασία επαναλαμβάνεται. Τα πεδία είναι χρονικά αποσυνδεδεμένα, έτσι ώστε να εξελίσσονται με ένα απλό, άμεσο σχήμα, όπου κάθε ενημέρωση εξαρτάται μόνο από τις τοπικές γειτονικές τιμές. Αυτό απεικονίζεται σχηματικά στο ?? και ουσιαστικά αυτό είναι ό,τι κάνει ο υπολογιστής σε κάθε χρονικό βήμα.

5. Ένα απλό παράδειγμα 1D FDTD σε Matlab/Python.

Το ίδιο μονοδιάστατο σχήμα FDTD μπορεί να υλοποιηθεί σε μόλις λίγες δεκάδες γραμμές κώδικα. Παρακάτω παρουσιάζονται ισοδύναμες υλοποιήσεις σε MATLAB και Python. Και οι δύο δημιουργούν έναν παλμό Gaussian και απεικονίζουν την ανάκλαση από μια διαφορά στον δείκτη διάθλασης.

clear;
clc;

% Physical constants
c0   = 299792458;
mu0  = 4*pi*1e-7;
eps0 = 1/(mu0*c0^2);

% Grid
Nx = 400;
dx = 20e-9;
dt = dx/(2*c0);

% Material profile
n = ones(1, Nx);
n(220:end) = 3.0;
eps_r = n.^2;

% Fields
Ez = zeros(1, Nx);
Hy = zeros(1, Nx-1);

% Source
source_pos = 80;
t0 = 40;
spread = 12;

Nt = 600;

figure;

for t = 1:Nt

    % Update H
    for i = 1:Nx-1
        Hy(i) = Hy(i) + (dt/(mu0*dx))*(Ez(i+1) - Ez(i));
    end

    % Update E
    for i = 2:Nx-1
        Ez(i) = Ez(i) + (dt/(eps0*eps_r(i)*dx))*(Hy(i) - Hy(i-1));
    end

    % Source
    Ez(source_pos) = Ez(source_pos) + exp(-((t - t0)/spread)^2);

    % Simple absorbing boundaries
    Ez(1)   = Ez(2);
    Ez(end) = Ez(end-1);

    % Plot every few steps
    if mod(t,5) == 0
        plot(1:Nx, Ez, 'b', 'LineWidth', 1.5);
        hold on;
        plot([220 220], [-1.5 1.5], '--r', 'LineWidth', 1.0);
        hold off;
        xlabel('Grid index');
        ylabel('E_z');
        title(['1D FDTD pulse reflecting from an index step, t = ', num2str(t)]);
        axis([1 Nx -1.5 1.5]);
        grid on;
        drawnow;
    end
end

import numpy as np
import matplotlib.pyplot as plt

c0 = 299792458.0
mu0 = 4.0e-7 * np.pi
eps0 = 1.0 / (mu0 * c0**2)

Nx = 400
dx = 20e-9
dt = dx / (2 * c0)

n = np.ones(Nx)
n[220:] = 3.0
eps_r = n**2

Ez = np.zeros(Nx)
Hy = np.zeros(Nx - 1)

source_pos = 80
t0 = 40
spread = 12

Nt = 600

plt.ion()
fig, ax = plt.subplots()

for t in range(Nt):

    Hy += (dt / (mu0 * dx)) * (Ez[1:] - Ez[:-1])
    Ez[1:-1] += (dt / (eps0 * eps_r[1:-1] * dx)) * (Hy[1:] - Hy[:-1])

    Ez[source_pos] += np.exp(-((t - t0) / spread) ** 2)

    Ez[0] = Ez[1]
    Ez[-1] = Ez[-2]

    if t % 5 == 0:
        ax.clear()
        ax.plot(Ez)
        ax.axvline(220, linestyle='--')
        ax.set_xlim(0, Nx)
        ax.set_ylim(-1.5, 1.5)
        plt.pause(0.01)

plt.ioff()
plt.show()

6. Εκτέλεση της προσομοίωσης.

Η συμπεριφορά της προσομοίωσης φαίνεται στα ??–??. Αρχικά, ο παλμός Gaussian διαδίδεται ομαλά μέσα στην περιοχή \(n=1\). Όταν φτάνει στο σημείο αλλαγής του δείκτη διάθλασης, το πεδίο διαχωρίζεται σε δύο συνιστώσες: έναν ανακλώμενο παλμό που ταξιδεύει πίσω στην περιοχή χαμηλού δείκτη διάθλασης και έναν διαθλώμενο παλμό που εισέρχεται στην περιοχή υψηλού δείκτη διάθλασης. Επειδή ο δείκτης διάθλασης αυξάνεται στο \(n=3\), το διαθλώμενο κύμα επιβραδύνεται και συμπιέζεται σε μήκος κύματος, ενώ η ανακλώμενη συνιστώσα διατηρεί τα αρχικά χαρακτηριστικά διάδοσης. Αυτό είναι ακριβώς η αναμενόμενη φυσική συμπεριφορά και αποδεικνύει ότι ακόμη και αυτή η ελάχιστη υλοποίηση αποτυπώνει την ουσία της φυσικής της διάδοσης, της ανάκλασης και της μετάδοσης των κυμάτων.

11. Πού να πάμε μετά;

Τα παραδείγματα MATLAB και Python που αναφέρθηκαν παραπάνω δείχνουν πόσο απλή είναι η αρχή της μεθόδου FDTD, αλλά πρόκειται μόνο για απλοποιημένες υλοποιήσεις. Στην πράξη, οι ρεαλιστικές οπτικές προσομοιώσεις απαιτούν μεγάλα πλέγματα, μεγάλους χρόνους εκτέλεσης και προσεκτικό χειρισμό των ορίων και των πηγών. Για να χρησιμοποιηθεί η μέθοδος FDTD σε χρήσιμες, πραγματικές εφαρμογές, χρειάζεστε έναν επιλυτή που να είναι ιδιαίτερα βελτιστοποιημένος, να υποστηρίζει παράλληλη επεξεργασία και να μπορεί να εκτελείται τόσο σε CPU όσο και σε GPU.

Αυτό ακριβώς προσφέρει η OghmaNano: μια ολοκληρωμένη τρισδιάστατη αριθμητική λύση FDTD, γραμμένη σε C με επιτάχυνση μέσω OpenCL, η οποία επιτρέπει στην ίδια ακριβώς μοντελοποίηση να εκτελείται απρόσκοπτα σε πολυεπεξεργαστικές CPU ή GPU. Αντί να δημιουργήσετε τη δική σας αριθμητική λύση από την αρχή, μπορείτε να χρησιμοποιήσετε τις ιδέες που παρουσιάζονται εδώ και να τις εφαρμόσετε άμεσα σε ρεαλιστικές φωτονικές δομές, συσκευές και προβλήματα κυματικής οπτικής.

Για να ξεκινήσετε, μπορείτε να εξερευνήσετε τη σελίδα FDTD στο OghmaNano για μια επισκόπηση των δυνατοτήτων, να επανεξετάσετε τι είναι η μέθοδος FDTD ή την πλήρη μαθηματική απόδειξη, ή να προχωρήσετε απευθείας στην πρακτική ρύθμιση χρησιμοποιώντας πηγές φωτός, οριακές συνθήκες και ανιχνευτές.

Επιπλέον, διατίθεται μια μεγάλη ποικιλία έτοιμων παραδειγμάτων, όπως διάδοση σε ελεύθερο χώρο, διπλή σχισμή και φαινόμενο συμβολής, φωτονικές κρυσταλλικές οπτικές ίνες, πλέγματα Bragg, κοιλότητες Fabry–Pérot, κεκλιμένες όψεις οπτικών ινών, κωνικές δομές, διατάξεις Mach–Zehnder, δακτυλιοειδείς συντονιστές και αντι-ανακλαστικές επιστρώσεις.

Εν ολίγοις: η βασική μέθοδος είναι απλή, όπως φαίνεται παραπάνω, αλλά το OghmaNano έχει σχεδιαστεί για να αξιοποιήσει αυτή την απλότητα και να την επεκτείνει σε ένα πρακτικό εργαλείο υψηλής απόδοσης για ρεαλιστικές οπτικές προσομοιώσεις.