Home Esempi Immagini dello schermo Manuale utente Logo di Bluesky YouTube ☰
OghmaNano Piattaforma di simulazione multifisica per dispositivi optoelettronici e sistemi fotonici SCARICA Guida rapida per iniziare.

Cos'è FDTD?Cos'è FDTD?

1. Introduzione1. Introduzione

FDTD simulation of a ring resonator coupled to a bus waveguide showing field intensity distribution.
Intensità del campo in un risonatore ad anello accoppiato a una guida d’onda lineare.Intensità del campo in un risonatore ad anello accoppiato a una guida d’onda lineare.
Three-dimensional FDTD simulation visualising energy density build-up in a resonant photonic structure.
Immagine tridimensionale ottenuta con il metodo FDTD che mostra l’accumulo di densità energetica in una struttura risonante.Immagine tridimensionale ottenuta con il metodo FDTD che mostra l’accumulo di densità energetica in una struttura risonante.

Il metodo delle differenze finite nel dominio del tempo (FDTD) risolve direttamente le equazioni di Maxwell – che descrivono la luce – nel dominio del tempo, utilizzando una griglia spaziale discreta. A differenza del ray tracing, dei metodi di risoluzione delle modalità o degli approcci basati sulla matrice di trasferimento, il FDTD non fa assunzioni semplificative sulla propagazione o sulla struttura modale. Invece, i campi elettrici e magnetici vengono calcolati esplicitamente nel tempo, discretizzando le equazioni del rotore – la legge di Ampère e la legge di Faraday – in modo che le onde elettromagnetiche si propaghino, si disperdano, interferiscano e si diffondano naturalmente in tutto il dominio di simulazione. Questa versatilità rende il FDTD un metodo molto preciso e flessibile, ma anche piuttosto impegnativo in termini di risorse di memoria e tempo di calcolo.Il metodo delle differenze finite nel dominio del tempo (FDTD) risolve direttamente le equazioni di Maxwell – che descrivono la luce – nel dominio del tempo, utilizzando una griglia spaziale discreta. A differenza del ray tracing, dei metodi di risoluzione delle modalità o degli approcci basati sulla matrice di trasferimento, il FDTD non fa assunzioni semplificative sulla propagazione o sulla struttura modale. Invece, i campi elettrici e magnetici vengono calcolati esplicitamente nel tempo, discretizzando le equazioni del rotore – la legge di Ampère e la legge di Faraday – in modo che le onde elettromagnetiche si propaghino, si disperdano, interferiscano e si diffondano naturalmente in tutto il dominio di simulazione. Questa versatilità rende il FDTD un metodo molto preciso e flessibile, ma anche piuttosto impegnativo in termini di risorse di memoria e tempo di calcolo.

\[ \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}. \] \[ \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}. \]

Poiché il metodo FDTD non si basa su ipotesi specifiche riguardanti la direzione di propagazione, la struttura del dispositivo o il comportamento del campo, una singola simulazione può analizzare la risposta su un ampio spettro di frequenze, la diffrazione e l’interferenza, gli effetti del campo vicino, la geometria su scala sub-lunghezza d’onda e i fenomeni elettromagnetici transitori, il tutto all’interno di un unico modello. Questa versatilità è ciò che rende il metodo FDTD così efficace, ma comporta anche un costo: la griglia spaziale deve risolvere la più piccola scala di lunghezza fisica, l’intervallo di tempo è strettamente vincolato dalla stabilità numerica e le simulazioni tridimensionali diventano rapidamente complesse in termini di memoria e tempo di calcolo, anche per geometrie che sembrano di dimensioni contenute.Poiché il metodo FDTD non si basa su ipotesi specifiche riguardanti la direzione di propagazione, la struttura del dispositivo o il comportamento del campo, una singola simulazione può analizzare la risposta su un ampio spettro di frequenze, la diffrazione e l’interferenza, gli effetti del campo vicino, la geometria su scala sub-lunghezza d’onda e i fenomeni elettromagnetici transitori, il tutto all’interno di un unico modello. Questa versatilità è ciò che rende il metodo FDTD così efficace, ma comporta anche un costo: la griglia spaziale deve risolvere la più piccola scala di lunghezza fisica, l’intervallo di tempo è strettamente vincolato dalla stabilità numerica e le simulazioni tridimensionali diventano rapidamente complesse in termini di memoria e tempo di calcolo, anche per geometrie che sembrano di dimensioni contenute.

2. Perché non dovreste utilizzare il metodo FDTD...2. Perché non dovreste utilizzare il metodo FDTD...

Nella maggior parte dei problemi ottici pratici, il metodo FDTD non dovrebbe essere la prima soluzione a cui ricorrere. È un metodo che richiede molte risorse computazionali e, per molti sistemi, esistono approcci molto più efficienti. La progettazione di sistemi di lenti, come fotocamere o telescopi, la modellazione di sistemi di imaging su larga scala o la simulazione del percorso della luce attraverso assemblaggi ottici macroscopici sono attività che si possono svolgere in modo molto più efficace utilizzando il metodo del tracciamento dei raggi. Allo stesso modo, le strutture planari multistrato, come le celle solari a film sottile, gli stack OLED e i rivestimenti ottici, possono essere analizzate in modo molto più efficiente utilizzando i metodi della matrice di trasferimento. In questi casi, l’utilizzo del metodo FDTD comporta un notevole aumento del carico computazionale senza fornire ulteriori informazioni fisiche.Nella maggior parte dei problemi ottici pratici, il metodo FDTD non dovrebbe essere la prima soluzione a cui ricorrere. È un metodo che richiede molte risorse computazionali e, per molti sistemi, esistono approcci molto più efficienti. La progettazione di sistemi di lenti, come fotocamere o telescopi, la modellazione di sistemi di imaging su larga scala o la simulazione del percorso della luce attraverso assemblaggi ottici macroscopici sono attività che si possono svolgere in modo molto più efficace utilizzando il metodo del tracciamento dei raggi. Allo stesso modo, le strutture planari multistrato, come le celle solari a film sottile, gli stack OLED e i rivestimenti ottici, possono essere analizzate in modo molto più efficiente utilizzando i metodi della matrice di trasferimento. In questi casi, l’utilizzo del metodo FDTD comporta un notevole aumento del carico computazionale senza fornire ulteriori informazioni fisiche.

La maggior parte dei metodi di simulazione elettromagnetica si basa, in ultima analisi, sulle equazioni di Maxwell. Tecniche come i metodi di risoluzione delle modalità e i metodi della matrice di trasferimento vengono sviluppate introducendo ipotesi semplificative sulla geometria, sulla dimensionalità o sulla struttura del campo, al fine di accelerare i calcoli. Queste approssimazioni sono estremamente efficaci, ma solo quando le ipotesi su cui si basano sono valide.La maggior parte dei metodi di simulazione elettromagnetica si basa, in ultima analisi, sulle equazioni di Maxwell. Tecniche come i metodi di risoluzione delle modalità e i metodi della matrice di trasferimento vengono sviluppate introducendo ipotesi semplificative sulla geometria, sulla dimensionalità o sulla struttura del campo, al fine di accelerare i calcoli. Queste approssimazioni sono estremamente efficaci, ma solo quando le ipotesi su cui si basano sono valide.

Il metodo FDTD si colloca all’estremo opposto di questo spettro. È estremamente efficace, in particolare per sistemi tridimensionali complessi, geometrie irregolari e problemi in cui il campo non può essere scomposto in modalità o livelli semplici, come ad esempio cristalli fotonici, strutture sub-lunghezza d’onda, mezzi che disperdono fortemente la luce e accoppiamenti di campo vicino. Il prezzo di questa flessibilità è rappresentato dai costi computazionali, sia in termini di memoria che di tempo di elaborazione.Il metodo FDTD si colloca all’estremo opposto di questo spettro. È estremamente efficace, in particolare per sistemi tridimensionali complessi, geometrie irregolari e problemi in cui il campo non può essere scomposto in modalità o livelli semplici, come ad esempio cristalli fotonici, strutture sub-lunghezza d’onda, mezzi che disperdono fortemente la luce e accoppiamenti di campo vicino. Il prezzo di questa flessibilità è rappresentato dai costi computazionali, sia in termini di memoria che di tempo di elaborazione.

L’aspetto fondamentale è scegliere lo strumento giusto. Utilizzate modelli semplificati quando le leggi della fisica lo consentono e utilizzate il metodo FDTD quando queste ipotesi non sono più valide. Il punto di forza di questo metodo è la sua flessibilità e la sua completezza fisica; il compromesso è rappresentato dal costo computazionale.L’aspetto fondamentale è scegliere lo strumento giusto. Utilizzate modelli semplificati quando le leggi della fisica lo consentono e utilizzate il metodo FDTD quando queste ipotesi non sono più valide. Il punto di forza di questo metodo è la sua flessibilità e la sua completezza fisica; il compromesso è rappresentato dal costo computazionale.

3. Le equazioni di Maxwell alla base del metodo FDTD3. Le equazioni di Maxwell alla base del metodo FDTD

In sostanza, il metodo FDTD è molto semplice. Si partono dalle equazioni di Maxwell, le si espandono nelle loro componenti per ciascuna direzione spaziale e quindi si risolvono direttamente su una griglia. Non c’è nulla di misterioso: il metodo è concettualmente semplice e, in molti aspetti, più facile da applicare rispetto ad approcci numerici più complessi, come i risolutori basati sui modelli di deriva-diffusione. Per iniziare, espandiamo la legge di Ampère. In un mezzo lineare con conducibilità \(\sigma\), permittività \(\epsilon\) e permeabilità \(\mu\), questa può essere scritta come:In sostanza, il metodo FDTD è molto semplice. Si partono dalle equazioni di Maxwell, le si espandono nelle loro componenti per ciascuna direzione spaziale e quindi si risolvono direttamente su una griglia. Non c’è nulla di misterioso: il metodo è concettualmente semplice e, in molti aspetti, più facile da applicare rispetto ad approcci numerici più complessi, come i risolutori basati sui modelli di deriva-diffusione. Per iniziare, espandiamo la legge di Ampère. In un mezzo lineare con conducibilità \(\sigma\), permittività \(\epsilon\) e permeabilità \(\mu\), questa può essere scritta come:

\[ \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 \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} \]

che può essere ampliato comeche può essere ampliato come

\[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z} \] \[ \sigma E_x + \epsilon \frac{\partial E_x}{\partial t} = \frac{\partial H_z}{\partial y} - \frac{\partial H_y}{\partial z} \]

\[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = -\frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z} \] \[ \sigma E_y + \epsilon \frac{\partial E_y}{\partial t} = -\frac{\partial H_z}{\partial x} + \frac{\partial H_x}{\partial z} \]

\[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} \] \[ \sigma E_z + \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y} \]

La legge di Faraday è espressa come segue:La legge di Faraday è espressa come segue:

\[ -\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} \] \[ -\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} \]

che può essere ampliato per ottenere:che può essere ampliato per ottenere:

\[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z} \] \[ -\sigma_m H_x - \mu \frac{\partial H_x}{\partial t} = \frac{\partial E_z}{\partial y} - \frac{\partial E_y}{\partial z} \]

\[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = -\frac{\partial E_z}{\partial x} + \frac{\partial E_x}{\partial z} \] \[ -\sigma_m H_y - \mu \frac{\partial H_y}{\partial t} = -\frac{\partial E_z}{\partial x} + \frac{\partial E_x}{\partial z} \]

\[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y} \] \[ -\sigma_m H_z - \mu \frac{\partial H_z}{\partial t} = \frac{\partial E_y}{\partial x} - \frac{\partial E_x}{\partial y} \]

Le equazioni evidenziate in blu e in rosso rappresentano la definizione del metodo FDTD e ne costituiscono il nucleo. Non esiste nulla di più complesso. Tutto il resto riguarda l’implementazione: la sostituzione delle derivate con differenze finite, la scelta delle condizioni al contorno, la definizione delle sorgenti luminose (eccitazioni), l’estrazione dei risultati mediante rivelatori e la scrittura di un codice efficiente (spesso parallelo) per far evolvere i campi nel tempo. Questi componenti aggiuntivi rendono più complesso un risolutore FDTD pratico, ma non alterano la semplicità intrinseca del metodo di base.Le equazioni evidenziate in blu e in rosso rappresentano la definizione del metodo FDTD e ne costituiscono il nucleo. Non esiste nulla di più complesso. Tutto il resto riguarda l’implementazione: la sostituzione delle derivate con differenze finite, la scelta delle condizioni al contorno, la definizione delle sorgenti luminose (eccitazioni), l’estrazione dei risultati mediante rivelatori e la scrittura di un codice efficiente (spesso parallelo) per far evolvere i campi nel tempo. Questi componenti aggiuntivi rendono più complesso un risolutore FDTD pratico, ma non alterano la semplicità intrinseca del metodo di base.

4. Un esempio di simulazione 1D FDTD applicata a un giocattolo.4. Un esempio di simulazione 1D FDTD applicata a un giocattolo.

Flowchart of the FDTD update loop showing field initialisation, magnetic-field update, electric-field update, boundary conditions, and output
Diagramma di flusso del ciclo di aggiornamento FDTD.Diagramma di flusso del ciclo di aggiornamento FDTD.

Per dimostrare quanto sia semplice il metodo FDTD, sia dal punto di vista concettuale che pratico, ridurremo ora le equazioni complete a un modello semplificato monodimensionale, implementabile in poche decine di righe di codice MATLAB o Python. L’obiettivo qui non è modificare i principi fisici, ma eliminare le complessità superflue, in modo che la struttura del metodo sia completamente trasparente.Per dimostrare quanto sia semplice il metodo FDTD, sia dal punto di vista concettuale che pratico, ridurremo ora le equazioni complete a un modello semplificato monodimensionale, implementabile in poche decine di righe di codice MATLAB o Python. L’obiettivo qui non è modificare i principi fisici, ma eliminare le complessità superflue, in modo che la struttura del metodo sia completamente trasparente.

Per fare ciò, consideriamo la variazione solo lungo la direzione \(x\) e impostiamo esplicitamente a zero le derivate spaziali rimanenti: \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. In questo modo, otteniamo un sistema elettromagnetico trasversale con campo elettrico \(E_z(x,t)\) e campo magnetico \(H_y(x,t)\).Per fare ciò, consideriamo la variazione solo lungo la direzione \(x\) e impostiamo esplicitamente a zero le derivate spaziali rimanenti: \[ \frac{\partial}{\partial y} = 0, \qquad \frac{\partial}{\partial z} = 0. \]. In questo modo, otteniamo un sistema elettromagnetico trasversale con campo elettrico \(E_z(x,t)\) e campo magnetico \(H_y(x,t)\).

Alla luce di queste ipotesi, le equazioni di Maxwell sulla divergenza si riducono a:Alla luce di queste ipotesi, le equazioni di Maxwell sulla divergenza si riducono a:

\[ \epsilon \frac{\partial E_z}{\partial t} = \frac{\partial H_y}{\partial x} \] \[ \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}. \] \[ \mu \frac{\partial H_y}{\partial t} = \frac{\partial E_z}{\partial x}. \]

Queste sono le equazioni continue fondamentali per il caso monodimensionale. Tuttavia, nella forma attuale, non possono essere risolte direttamente su un computer. Per farlo, sostituiamo le derivate con differenze finite su una griglia discreta. Questo è il principio fondamentale del metodo: per una panoramica generale, consultare differenze finite.Queste sono le equazioni continue fondamentali per il caso monodimensionale. Tuttavia, nella forma attuale, non possono essere risolte direttamente su un computer. Per farlo, sostituiamo le derivate con differenze finite su una griglia discreta. Questo è il principio fondamentale del metodo: per una panoramica generale, consultare differenze finite.

Introducendo una griglia spaziale \(x_i = i\Delta x\) e intervalli temporali \(t^n = n\Delta t\) e utilizzando la disposizione standard a scacchiera (leapfrog), le equazioni di aggiornamento diventano:Introducendo una griglia spaziale \(x_i = i\Delta x\) e intervalli temporali \(t^n = n\Delta t\) e utilizzando la disposizione standard a scacchiera (leapfrog), le equazioni di aggiornamento diventano:

\[ 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] \] \[ 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]. \] \[ 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]. \]

Questa è già la versione completa del metodo FDTD nella sua forma più semplice. Il campo elettrico viene aggiornato in base alle differenze spaziali del campo magnetico, e il campo magnetico viene aggiornato in base alle differenze spaziali del campo elettrico, con i due campi che si alternano nel tempo.Questa è già la versione completa del metodo FDTD nella sua forma più semplice. Il campo elettrico viene aggiornato in base alle differenze spaziali del campo magnetico, e il campo magnetico viene aggiornato in base alle differenze spaziali del campo elettrico, con i due campi che si alternano nel tempo.

Per un materiale non magnetico, \(\mu = \mu_0\), la permittività è correlata all’indice di rifrazione tramite la seguente relazione:Per un materiale non magnetico, \(\mu = \mu_0\), la permittività è correlata all’indice di rifrazione tramite la seguente relazione:

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

Questo significa che una variazione dell’indice di rifrazione viene introdotta semplicemente modificando il valore di \(\epsilon\) in diversi punti della griglia. Nell’esempio riportato di seguito, consideriamo un mezzo che passa da \(n=1\) a \(n=3\), in modo che un impulso emesso dalla parte sinistra venga parzialmente riflesso e parzialmente trasmesso all’interfaccia.Questo significa che una variazione dell’indice di rifrazione viene introdotta semplicemente modificando il valore di \(\epsilon\) in diversi punti della griglia. Nell’esempio riportato di seguito, consideriamo un mezzo che passa da \(n=1\) a \(n=3\), in modo che un impulso emesso dalla parte sinistra venga parzialmente riflesso e parzialmente trasmesso all’interfaccia.

In pratica, le equazioni di cui sopra vengono semplicemente inserite all’interno di un ciclo iterativo che simula l’evoluzione nel tempo: il campo magnetico viene aggiornato a partire dal campo elettrico, vengono applicate le condizioni al contorno, quindi il campo elettrico viene aggiornato a partire dal campo magnetico, e il processo si ripete. I campi sono sfalsati nel tempo in modo da avanzare in modo semplice e diretto, secondo uno schema esplicito, e ogni aggiornamento dipende esclusivamente dai valori locali dei punti adiacenti. Questo è illustrato schematicamente in ?? e rappresenta essenzialmente tutto ciò che il computer esegue a ogni passo temporale.In pratica, le equazioni di cui sopra vengono semplicemente inserite all’interno di un ciclo iterativo che simula l’evoluzione nel tempo: il campo magnetico viene aggiornato a partire dal campo elettrico, vengono applicate le condizioni al contorno, quindi il campo elettrico viene aggiornato a partire dal campo magnetico, e il processo si ripete. I campi sono sfalsati nel tempo in modo da avanzare in modo semplice e diretto, secondo uno schema esplicito, e ogni aggiornamento dipende esclusivamente dai valori locali dei punti adiacenti. Questo è illustrato schematicamente in ?? e rappresenta essenzialmente tutto ciò che il computer esegue a ogni passo temporale.

5. Un esempio di simulazione 1D FDTD per giocattoli, realizzato in Matlab/Python.5. Un esempio di simulazione 1D FDTD per giocattoli, realizzato in Matlab/Python.

Lo stesso schema FDTD monodimensionale può essere implementato con poche decine di righe di codice. Di seguito sono riportate implementazioni equivalenti in MATLAB e Python. Entrambe simulano l’emissione di un impulso gaussiano e mostrano la riflessione che si verifica in corrispondenza di una discontinuità dell’indice di rifrazione.Lo stesso schema FDTD monodimensionale può essere implementato con poche decine di righe di codice. Di seguito sono riportate implementazioni equivalenti in MATLAB e Python. Entrambe simulano l’emissione di un impulso gaussiano e mostrano la riflessione che si verifica in corrispondenza di una discontinuità dell’indice di rifrazione.

MATLABMATLAB


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

Python (NumPy/Matplotlib)Python (NumPy/Matplotlib)


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. Esecuzione della simulazione6. Esecuzione della simulazione

1D FDTD pulse at early time showing initial propagation before reaching the index step
Propagazione nelle prime fasi (\(t \approx 110\)).Propagazione nelle prime fasi (\(t \approx 110\)).
1D FDTD pulse interacting with an index step showing partial reflection and transmission
Interazione con la fase di indicizzazione (\(t \approx 280\)).Interazione con la fase di indicizzazione (\(t \approx 280\)).
1D FDTD pulse after interaction with index step showing reflected and transmitted waves
Impulsi riflessi e trasmessi (\(t \approx 570\)).Impulsi riflessi e trasmessi (\(t \approx 570\)).

Il comportamento della simulazione è illustrato nelle figure ??–??. Inizialmente, l’impulso gaussiano si propaga in modo uniforme attraverso la regione \(n=1\). Quando raggiunge il punto di discontinuità dell’indice di rifrazione, il campo si divide in due componenti: un impulso riflesso che si propaga all’indietro nella regione a basso indice e un impulso trasmesso che entra nella regione a indice più elevato. Poiché l’indice di rifrazione aumenta fino a \(n=3\), l’onda trasmessa rallenta e si comprime in lunghezza d’onda, mentre la componente riflessa mantiene le caratteristiche di propagazione originali. Questo corrisponde esattamente al comportamento fisico previsto e dimostra che anche questa implementazione minima riproduce la fisica essenziale della propagazione, della riflessione e della trasmissione delle onde.Il comportamento della simulazione è illustrato nelle figure ??–??. Inizialmente, l’impulso gaussiano si propaga in modo uniforme attraverso la regione \(n=1\). Quando raggiunge il punto di discontinuità dell’indice di rifrazione, il campo si divide in due componenti: un impulso riflesso che si propaga all’indietro nella regione a basso indice e un impulso trasmesso che entra nella regione a indice più elevato. Poiché l’indice di rifrazione aumenta fino a \(n=3\), l’onda trasmessa rallenta e si comprime in lunghezza d’onda, mentre la componente riflessa mantiene le caratteristiche di propagazione originali. Questo corrisponde esattamente al comportamento fisico previsto e dimostra che anche questa implementazione minima riproduce la fisica essenziale della propagazione, della riflessione e della trasmissione delle onde.

11. Qual è la prossima destinazione?11. Qual è la prossima destinazione?

Gli esempi MATLAB e Python riportati sopra illustrano quanto sia semplice, in linea di principio, il metodo FDTD, ma si tratta solo di implementazioni di base. Nella pratica, simulazioni ottiche realistiche richiedono griglie ampie, tempi di calcolo prolungati e una gestione accurata dei bordi e delle sorgenti. Per applicare il metodo FDTD a problemi reali e concreti, è necessario un risolutore altamente ottimizzato, parallelizzato e in grado di funzionare sia su CPU che su GPU.Gli esempi MATLAB e Python riportati sopra illustrano quanto sia semplice, in linea di principio, il metodo FDTD, ma si tratta solo di implementazioni di base. Nella pratica, simulazioni ottiche realistiche richiedono griglie ampie, tempi di calcolo prolungati e una gestione accurata dei bordi e delle sorgenti. Per applicare il metodo FDTD a problemi reali e concreti, è necessario un risolutore altamente ottimizzato, parallelizzato e in grado di funzionare sia su CPU che su GPU.

OghmaNano offre proprio questo: un potente risolutore 3D FDTD, scritto in linguaggio C e dotato di accelerazione OpenCL, che consente di eseguire lo stesso modello in modo ottimale su CPU o GPU multi-core. Invece di sviluppare un risolutore personalizzato da zero, è possibile utilizzare le idee presentate qui e applicarle immediatamente a strutture fotoniche, dispositivi e problemi di ottica ondulatoria realistici.OghmaNano offre proprio questo: un potente risolutore 3D FDTD, scritto in linguaggio C e dotato di accelerazione OpenCL, che consente di eseguire lo stesso modello in modo ottimale su CPU o GPU multi-core. Invece di sviluppare un risolutore personalizzato da zero, è possibile utilizzare le idee presentate qui e applicarle immediatamente a strutture fotoniche, dispositivi e problemi di ottica ondulatoria realistici.

Per iniziare, potete consultare la pagina FDTD in OghmaNano per avere una panoramica delle funzionalità, oppure ripassare cos'è FDTD o la derivazione matematica completa, oppure passare direttamente alla configurazione pratica utilizzando sorgenti di luce, condizioni al contorno e rivelatori.Per iniziare, potete consultare la pagina FDTD in OghmaNano per avere una panoramica delle funzionalità, oppure ripassare cos'è FDTD o la derivazione matematica completa, oppure passare direttamente alla configurazione pratica utilizzando sorgenti di luce, condizioni al contorno e rivelatori.

È disponibile anche un’ampia gamma di esempi già pronti per l’uso, tra cui: propagazione in spazio libero, interferenza a doppia fenditura, guide d’onda a cristallo fotonico, reticoli di Bragg, cavità di Fabry–Pérot, facce inclinate delle guide d’onda, strutture a forma di cono, dispositivi Mach–Zehnder, risuonatori ad anello e rivestimenti antiriflesso.È disponibile anche un’ampia gamma di esempi già pronti per l’uso, tra cui: propagazione in spazio libero, interferenza a doppia fenditura, guide d’onda a cristallo fotonico, reticoli di Bragg, cavità di Fabry–Pérot, facce inclinate delle guide d’onda, strutture a forma di cono, dispositivi Mach–Zehnder, risuonatori ad anello e rivestimenti antiriflesso.

In sintesi: il metodo di base è semplice, come illustrato in precedenza, ma OghmaNano è progettato per sfruttare questa semplicità e ampliarla, trasformandola in uno strumento pratico e ad alte prestazioni per la simulazione ottica reale.In sintesi: il metodo di base è semplice, come illustrato in precedenza, ma OghmaNano è progettato per sfruttare questa semplicità e ampliarla, trasformandola in uno strumento pratico e ad alte prestazioni per la simulazione ottica reale.