Un milione di elementi, un minuto.
Abbiamo risolto un problema di correnti parassite (eddy currents) con un milione di elementi conduttori su un laptop con 16 GB di RAM, alla frequenza di 1 MHz, in circa un minuto, utilizzando circa 4 GB di memoria. Un milione di elementi rappresenta esclusivamente il rame: ogni pista e ogni spira di filo. L'aria circostante non viene mai discretizzata, ma il campo elettromagnetico nello spazio circostante può comunque essere ricavato dalla soluzione.
A 1 MHz, la profondità di penetrazione (skin depth) nel rame è di circa 65 micrometri, quasi il doppio dei 35 micrometri di spessore del rame da 1 oz. A questa frequenza, quindi, la corrente sta già iniziando a ridistribuirsi all'interno della sezione trasversale della pista.Il costo computazionale della soluzione non dipende dalla frequenza, nell'intero intervallo in cui è valida l'approssimazione magnetoquasistatica. Al di sopra di tale intervallo è necessario adottare un modello elettromagnetico a onda completa (full-wave): si tratta di un limite del modello fisico, non del solutore.
Elementi finiti
Ho vent'anni di esperienza nella risoluzione di problemi di correnti parassite con il metodo degli elementi finiti, compreso lo sviluppo dei relativi codici di calcolo.
Non si tratta quindi di una critica al software di altri. Tre aspetti impediscono a questo metodo di affrontare efficacemente problemi di questa classe.
Discretizzazione dell'aria.
Anche l'aria deve essere discretizzata. Nel caso di un circuito stampato, un milione di elementi nel rame significa da dieci a venti milioni di elementi nell'aria, spesso anche di più. In molti casi reali, costruire la mesh dell'aria attorno a una bobina multistrato con rame spesso 35 micrometri è semplicemente impossibile.
Convergenza.
Nell'aria la conducibilità è nulla, quindi scompare il termine che rende non singolare il sistema curl-curl, lasciando un nucleo di gradienti di dimensioni molto elevate. La gravità del problema dipende dalla formulazione e dal precondizionatore utilizzati. Nelle geometrie dominate dall'aria non abbiamo trovato alcuna combinazione, nemmeno con metodi multigrid, capace di convergere in poche iterazioni.
Movimento.
Generalmente il movimento richiede di ricostruire la mesh per ogni posizione, a meno che non presenti una simmetria sfruttabile tramite un'interfaccia scorrevole (sliding interface). Questa soluzione funziona per una macchina rotante, ma non per un target che, oltre a ruotare, si inclina e si sposta.
Metodi integrali e cosa comprende questa definizioneI metodi integrali si suddividono in tre categorie: metodi integrali di volume, metodi integrali di superficie e PEEC (Partial Element Equivalent Circuit).
Sviluppato utilizzando parallelepipedi con direzioni della corrente fissate a priori. Introducendo funzioni di base analoghe a quelle degli elementi finiti, diventa un metodo integrale di volume. Tutti condividono la proprietà fondamentale per il problema che stiamo affrontando: vengono discretizzati soltanto i conduttori, senza mesh dell'aria né superfici artificiali di troncamento del dominio. Ciò che li distingue è il modo in cui gestiscono la matrice densa che deriva da questa scelta.
Matrici dense
La matrice è piena, quindi la memoria necessaria cresce con il quadrato del numero di incognite. Nel benchmark con un milione di elementi pubblicato nel nostro articolo su CMAME, la formulazione standard con matrice densa richiede 8.141 GB, circa otto terabyte.Non c'è alcun parametro da ottimizzare che possa risolvere il problema.
ACA – Adaptive Cross Approximation
È la risposta standard proposta in letteratura, ma vale la pena essere espliciti su quanto siano limitati i benefici. Su un modello di 252.000 tetraedri, la compressione ha ridotto il picco di memoria all'8% rispetto al caso con matrice densa, passando da 426 GB a 34,8 GB. Il tempo totale, invece, è sceso soltanto al 43%, da 26.495 a 11.437 secondi. Un miglioramento di un fattore dodici nella memoria, ma soltanto di un fattore due nel tempo, per un problema in cui il divario da colmare è di tre ordini di grandezza. Aumentando la dimensione del problema, questo divario diventa decisivo. Nel benchmark con un milione di elementi, sulla stessa macchina, ACA ha impiegato oltre 24 ore, mentre il metodo descritto più avanti ha richiesto 649 secondi. Inoltre, ACA necessitava di più memoria anche per un problema di dimensioni pari a un quarto. La compressione permette di far entrare questi modelli in memoria. Non è ciò che li rende veloci.
Fast Multipole Method
Tra le tecniche di compressione pubblicate, è quella che ha ottenuto i risultati migliori nei nostri benchmark. L'obiezione abituale riguarda il numero di iterazioni, ma esiste un secondo problema che quasi mai viene considerato: il costo di una singola interazione. Con le usuali funzioni di base definite sugli elementi, l'interazione tra una coppia di elementi non richiede il calcolo di un solo integrale, ma di un integrale per ogni coppia di funzioni di base. Questo significa 36 integrali di interazione pesati per una base a sei spigoli e 16 per una base di Raviart–Thomas o Rao–Wilton–Glisson associata alle facce. Il framework VINCO utilizza funzioni di base uniformi su ciascun elemento, per cui l'interazione tra due tetraedri richiede un solo integrale. Non si tratta di sacrificare la precisione in cambio della velocità. Nei nostri articoli dimostriamo l'equivalenza delle due formulazioni: la discretizzazione VINCO produce la stessa soluzione discreta e, su una mesh tetraedrica, la sua matrice di rigidezza coincide con quella costruita utilizzando funzioni di base di Raviart–Thomas.Per essere precisi, tale fattore riguarda il costo della parte diretta dell'operatore, ossia le interazioni di campo vicino (near field). L'accelerazione Fast Multipole modifica invece il modo in cui vengono valutate le interazioni di campo lontano (far field).Non elimina la complessità introdotta dalle funzioni di base che variano nello spazio, ma il costo aggiuntivo residuo dipende dall'implementazione. I fattori 16 e 36 vanno quindi interpretati come costi relativi delle interazioni di campo vicino, non come fattori universali applicabili al costo di ogni iterazione.
FFT su voxel
Tra tutti gli approcci descritti, questo è l'unico realmente competitivo in termini di fattibilità computazionale.Il confronto, però, deve essere effettuato a parità di accuratezza. Ed è proprio qui che emergono i problemi. Un milione di voxel non equivale a un milione di tetraedri. Disponiamo di dati relativi a un caso reale: la bobina trasmittente di un sensore induttivo di posizione, con rame spesso 35 micrometri su un circuito stampato da 0,8 mm. Il dispositivo è stato misurato con un Bode100 nell'intervallo da 10 kHz a 10 MHz e simulato con tutti e tre gli approcci. Il metodo FFT richiedeva una griglia computazionale di circa 120 milioni di celle, con risoluzione di 35 micrometri, mentre la formulazione con mesh conforme alla geometria utilizzava circa 800.000 tetraedri conduttori. Il modello tetraedrico mostrava inoltre un accordo migliore con le misure sperimentali, proprio perché rappresentava correttamente le parti curve della geometria. C'è poi un secondo costo nascosto in quella griglia. La FFT opera sull'intero volume computazionale, non soltanto sulle celle contenenti rame. In quell'esempio, 229.000 celle conduttrici erano contenute in una griglia di 120 milioni di celle, con un rapporto tra volume d'aria e volume conduttore pari a 522.Il costo dell'aria rientra così dalla finestra, proprio in una famiglia di metodi nata per evitarlo.
Non si tratta di una critica costruita per un articolo divulgativo.
Il nostro articolo su CMAME individua tre limitazioni dell'approccio basato sui voxel:
- La FFT ottiene la propria velocità imponendo una mesh cartesiana.
- La geometria voxelizzata introduce il noto errore a gradini (staircase error), lo stesso che porta generalmente a preferire gli elementi finiti alle differenze finite.
- Una griglia di voxel non consente un raffinamento locale della mesh.La conclusione di quello studio, riportata fedelmente, è che le griglie di voxel accelerate mediante FFT sono adatte a problemi con meno di circa un milione di elementi, senza geometrie curve e con geometrie ben allineate alla griglia strutturata.
Una bobina su circuito stampato non soddisfa nessuna di queste tre condizioni.
