Tecniche di Precondizionamento Inverso in un Modello Oceanico
Creators: name: Farina, Raffaele
Publisher: name: Università degli Studi di Napoli Federico II
Year: 2011

description: Molte importanti applicazioni di interesse scientifico e industriale sono sempre
riconducibili a sistemi dinamici basati sulle equazioni differenziali alle
derivate parziali, tra queste ritroviamo anche i modelli oceanici globali.
I modelli oceanici sono una componente dei modelli climatici e sono in crescente
sviluppo essendo utilizzati, non solo per lo studio delle dinamiche oceaniche, ma per l'intero sistema climatico globale.
I modelli in questione forniscono strumenti per l'interpretazione delle osservazioni oceaniche, la valutazione di scenari futuri come ad esempio quelli legati al riscaldamento climatico indotto dall'uomo e infine strumenti per le previsioni delle dinamiche degli oceani su differenti scale spaziali e temporali.
I modelli oceanici si basano sulle equazioni primitive della fluidodinamica e le variabili di interesse sono il campo delle velocità, la temperatura, la salinità, la pressione e la densità del fluido oceanico. Le principali differenze tra i vari modelli attualmente in uso sono nella parametrizzazione di fenomeni fisici di piccola scala, nel calcolo del gradiente di pressione superficiale, nell’adozione del sistema di coordinate verticali, nella scelta delle griglie orizzontali per la rappresentazione delle quantità fisiche del modello e infine nell'uso di diversi metodi numerici per la discretizzazione delle equazioni evolutive nelle equazioni fondamentali della fluidodinamica.
Nei modelli oceanici le scale temporali di interesse variano da mesi ad anni, ma nelle simulazioni alcuni dei fenomeni richiedono scale temporali che vanno da secondi a ore. Inoltre le simulazioni dei modelli vengono eseguite sull'intero pianeta e non si conosce se vi siano luoghi sulla Terra dove la risoluzione delle griglie è meno importante. Conseguentemente i costi computazionali
di questi modelli sono enormemente grandi tanto da limitare l'uso a volte di alcuni di questi.
Il lavoro di tesi si inserisce in questo contesto, con l’obiettivo di affrontare nello specifico le problematiche numeriche relative ad un solutore iterativo, il
Gradiente Coniugato (CG), utilizzato per la risoluzione di un sistema lineare, derivante dal nucleo ellittico di un modello per la circolazione oceanica globale
chiamato NEMO (Nucleus for European Modeling of the Ocean) introdotto per "chiudere" il modello alle equazioni primitive. Un simile sistema ovvero un sistema avente sempre la stessa matrice incompleta è risolto ad ogni passo
temporale del modello oceanico. Sperimentalmente si è osservato una lenta convergenza del solutore iterativo
CG nel modello quando le risoluzioni delle griglie di discretizzazioni aumentano o quando il dominio di simulazione include luoghi in prossimità dei poli
terrestri. 
Il sistema lineare, Ax = b, che si vuole esaminare è ottenuto dal nucleo ellittico del modello oceanico globale NEMO per mezzo delle differenze finite nello spazio. La sua matrice dei coefficienti A gode della proprietà di
sparsità, di simmetria e definita positività. Per avere un' idea delle dimensioni di questo sistema si fa riferimento ad ORCA025, che è uno delle
principali configurazioni del modello, in cui l'ordine di grandezza per la matrice incompleta A è di un milione.
La presente tesi evidenzia in principio le cause principali della scarsa rapidità di convergenza del solutore dovute, in primo luogo, all'aumento della dimensione del sistema causato dall' incremento della risoluzione delle griglie
spaziali. Un' analisi più approfondita, invece, rivela che questa insufficienza della velocità di convergenza è provocata anche dalla variazione del rapporto delle funzioni α e β rappresentanti i coefficienti del problema ellittico in esame.
E' noto che, l'uso dei precondizionatori nei metodi iterativi si rivela essere fondamentale per aumentare la rapidità di convergenza di questi nella risoluzione di un sistema lineare Ax = b, con delle definite proprietà.
L'obiettivo è ottenere un'approssimazione della soluzione, nei noti spazi di Krylov, nel minor numero di iterazioni possibile al fine di ridurre il tempo
di risoluzione e gli effetti della aritmetica a precisione finita quando questi vengono eseguiti su calcolatore. Il precondizionatore utilizzato nel modello
oceanico NEMO è quello diagonale.
All'interno della tesi è presentato un primo risultato che lega la distanza della matrice precondizionata del sistema alla dimensione di A e al rapporto tra le funzioni α e β nella norma di Frobenius mediante una stima di un limite inferiore ed uno superiore per la grandezza in disamina. Questo risultato caratterizza appunto la lenta convergenza del CG quando la dimensione del sistema lineare aumenta
e quando il modello è risolto in prossimità dei poli terrestri a causa di uno sbilanciamento tra le funzioni α e β. Essendo A una matrice simmetrica e definita positiva, si introduce nella tesi il precondizionatore ottenuto dalla fattorizzazione incompleta di Cholesky. La decomposizione incompleta di Cholesky è un potente precondizionatore per problemi simmetrici, definiti positivi e con elevata sparsità. Si dimostra nel nostro caso, attraverso 
esperimenti numerici, come questo precondizionatore sia in grado di annullare gli effetti della aumentata dimensione del sistema lineare Ax = b, a seguito dell'intensificazione delle griglie di discretizzazione e dello sbilanciamento
delle funzioni α e β nelle adiacenze dei poli. Questo precondizionatore ha la facoltà di incrementare notevolmente la rapidità di convergenza del solutore
iterativo nel sistema lineare Ax = b del modello oceanico NEMO. Purtroppo per grandi dimensioni di A, lo svantaggio principale dell’uso del precondizionatore di Cholesky sta nella risoluzione di un sistema lineare, di matrice incompleta il precondizioantore di Cholesky, eseguita ad ogni iterazione del solutore iterativo ed alla sua difficoltà di parallelizzazione.
Per superare questa problematica lo studio si è rivolto ad un' altra classe di precondizionatori, chiamati precondizionatori inversi sparsi, i quali cercano
di approssimare in modo diretto l'inversa di A. La loro implementazione richiede all'interno del solutore l'impiego dell'algoritmo del matrice per vettore sparso per ogni iterazione del solutore.
Su questa nuova classe di precondizionatori, grazie a delle assunzioni sulla matrice A ed a delle stime pervenute sugli elementi dell'inversa del fattore di Cholesky  di A, si è costruito un nuovo precondizionatore inverso ad hoc
per il problema dato, capace di ridurre il tempo di risoluzione del PCG con precondizionatore diagonale nel caso dello sbilanciamento dei coefficienti α e
β del 25%.
Lo studio di opportuni precondizionatori che accelerano la rapidità  di convergenza dei solutori iterativi si inserisce, in modo significativo, anche nel contesto futuro che vede le comunità climatiche utilizzare le nuove architetture di computing, formate da sistemi ibridi con GPU (Graphical Processing Unit) a core paralleli e CPU (Control Processor Unit) multi-core. L' utilizzo
delle GPU NVIDIA, consente di trarre, in maniera economica, un enorme vantaggio grazie alla loro elevata capacità di elaborazione parallela, mediante l'utilizzo di centinaia di core, completamente programmabili e in grado di
eseguire migliaia di thread congiuntamente. Nella tesi, viene anche presentata un'applicazione per GPU NVIDIA, sviluppata in ambiente CUDA (Compute Unified Device Architecture) e CUBLAS (CUDA Basic Linear Algebra),
che implementa il metodo del PCG (Gradiente Coniugato Precondizionato) in singola precisione con precondizionatore diagonale P, per la risoluzione del nucleo ellittico in disamina. In termini di efficienza, per gli ordini di grandezza della dimensione del problema sperimentati, fissato il valore della tolleranza sul residuo alla singola precisione, la versione parallela del PCG ottenuta mediante CUDA e CUBLAS si è dimostrata capace di ridurre i tempi di elaborazione dell'algoritmo sequenziale fino a 6 volte con l' hardware utilizzato.
Infine il motivo principale della proposta dell'utilizzo di precondizioatori di tipo inversi sparsi per il sistema lineare del modello è dovuto non solo al fatto che hanno una elevata capacità di ridurre il numero di iterazioni del
solutore per raggiungere la convergeza, ma anche perchè la loro implementazione, richiede il calcolo del matrice per vettore ad ogni iterazione del solutore.
Questa operazione è caratterizzata da un elevato grado di parallelismo sui dati che la rende particolarmente adatto all' elaborazione parallela mediante
CUDA. L'uso dei precondizionatori inversi sparsi nel modello ha ridotto il tempo di risoluzione del PCG con precondizioantore diagonale su GPU fino all’ 87% nel caso dello sbilanciamento dei coefficienti α e β del nucleo
ellittico.
I risultati pervenuti in questa tesi vogliono essere una proposta all'uso di tecniche di precondizionamento più efficienti nei modelli oceanici, al fine di superare i problemi di scarsa convergenza dei metodi iterativi all'interno dei modelli oceanici globali. Inoltre, anche in vista del futuro uso, da parte della comunità scientifica oceanografica, delle rivoluzionarie architetture di calcolo
parallele GPU NVIDIA, i risultati in questo lavoro anticipano l'uso di precondizionatori inversi sparsi al fine di sfruttare al meglio il tipo di parallelismo
di tali architetture di calcolo.
Si termina ribadendo che accelerare i solutori iterativi nei modelli oceanici globali è un obiettivo fondamentale nel quadro futuro della modellazione oceanica per spingere oltre nel tempo le simulazioni per le previsioni Climatiche.