Tecniche di Precondizionamento Inverso in un Modello Oceanico

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.