Applicazioni di Calcolo Parallelo in Astrofisica

1. Applicazioni di Calcolo Parallelo in Astrofisica Claudio Gheller CINECA c.gheller@cineca.it

2. Astronomia e Computer: un legame molto stretto. DALLE OSSERVAZIONI ALLA TEORIA L’uso del calcolatore è indispensabile in ogni passo della ricerca astronomica: • Sviluppo e progettazione di nuovi strumenti di osservazione • Controllo e gestione delle strumentazioni (telescopi, satelliti) • Acquisizione dei dati, telemetria • Data management e compressione • Riduzione e analisi dei dati • Elaborazione dei risultati • Realizzazione di modelli teorici • Confronto tra osservazioni e modelli teorici

3. Simulazioni Numeriche Le simulazioni numeriche, sono la realizzazione dei modelli teorici • Descrizione del modello • Condizioni iniziali • Equazioni evolutive • Rappresentazione dello spazio • Approssimazione del modello • Algoritmi di soluzione delle equazioni • Ricostruzione degli osservabili • Confronto con le osservazioni

4. Simulazioni numeriche in astrofisica I campi applicativi sono innumerevoli: Il sistema solare Il sole, nascita ed evoluzione dei pianeti… Le stelle Struttura, evoluzione, ammassi… Le galassie Nascita, evoluzione… La cosmologia Proprietà ed evoluzione dell’Universo

5. Cosmological Simulations Cosmological simulations are challenging applications since: What is a Cosmological Simulation ??? A cosmological simulation has the object of describing the birth and the evolution of the cosmological structures (galaxies, galaxy clusters, filaments) following the dynamics of all the components: baryons (typically fluid description), dark matter and stars (N-body)… Lo scambio di dati è gestito da un unico processo, che specifica sia l’origine che la destinazione dei dati stessi (one-side communication)

6. Applicazioni Cosmologiche Viene simulato un volume di dimensioni confrontabili all’intero universo (1/100) Cosa possiamo imparare ??? Fisica fondamentale e universo primordiale Interazioni fondamentali, Big Bang, Inflazione, transizioni di fase… COBE mission: map of the infrared sky La più recente immagine del Big Bang (MAP)

7. Struttura e proprietà delle galassie e degli ammassi di galassie Cosa possiamo imparare??? X-ray imageof the Virgo Cluster 16 Mpc* faraway 2.5° radius (compare to moon 0,25°) X-ray imageof the Coma Cluster 100 Mpc faraway More than 1000 galaxies X-ray imageof a Simulated Cluster Standard CDM model 2563 simulation *1 Mpc = 1019 km Temperature of infalling matter Temperatures range from 0 to 108 K …e così via…

8. Quanto costa una simulazione ??? • le strutture cosmologiche sono influenzate da processi che si sviluppano su un range di scale enorme; alcuni esempi: • le stelle si formano su scale ~ 1 parsec (1016km) • le galassie sono su scale ~10 kparsec • le forze di pressione influenzano scale ~ 10 Mparsec • la gravità agisce sino a scale comparabili con quella dell'universo, ~ 1000 Mparsec • Il range dinamico ideale sarebbe quindi di 109 !!!!!!!! • Vediamo cosa significa questo in termini di memoria del calcolatore:

9. Descrizione del problema Le osservazioni ci dicono che l'universo e' costituito da 3 componenti fondamentali: la materia barionica: protoni, neutroni- costituisce tutto cio' che emette "luce" ; la radiazione: prodotta dai corpi celesti e residuo del Big Bang (radiazione di fondo) - la sua densita' di energia e' molto minore di quella della materia; la materia oscura: natura ignota (neutrini massivi, altre particelle elmentari, pianeti ???...) - scoperta nelle galassie confrontando le curve di rotazione delle stelle (indice della distribuzione di massa totale) con le curve di luminosita' (indice della distribuzione di barioni). Confermata poi da innumerevoli altre osservazioni

10. Soluzione del Problema • Focalizziamoci sulle due componenti materiali: • la materia oscura ha solo interazioni a lungo raggio (gravita'). Per essa e' adeguato un approccio numerico tipo N-Body: • la materia e' rappresentata come un insieme di particelle il cui comportamento e' descritto dalla soluzione delle relative equazioni del moto • la materia barionica risente anche di interazioni a corto raggio: • essa quindi puo' venire piu' adeguatamente descritta come un fluido,la cui evoluzione viene caratterizzata dalla soluzione di un sistema di equazioni che stabiliscono la conservazione di • materia • quantita' di moto • energia • Infine e' necessario calcolare il campo gravitazionale che accoppia le due componenti.

11. Equazioni del modello Materia oscura: equazioni del moto. Materia Barionica: equazioni di conservazione. + equazione di stato Campo e Forze Gravitazionali: equazione di Poisson.

12. L’algoritmo Fluidodinamico alle Differenze Finite • Con l'approccio alle differenze finite spazio e tempo vengono rappresentati in modo discreto, come griglia computazionale multidimensionale. • Su questa griglia le vengono definite le quantità caratteristiche del sistema. Queste possono venire interpretate come: • valori medi sulle cella della griglia, • valori puntuali sui nodi della griglia. • Sulla griglia spazio-temporale, vengono riscritte in modo approssimato le equazioni che caratterizzano la dinamica del sistema, in modo che esse dipendano dai valori delle quantità caratteristiche approssimati sulla griglia stessa.

13. L’algoritmo Fluidodinamico alle Differenze Finite • Con l'approccio alle differenze finite spazio e tempo vengono rappresentati in modo discreto, come griglia computazionale multidimensionale. • L'idea base del modello numerico e' la seguente: • La materia barionica viene descritta come un fluido; • il fluido e' descritto attraverso la sua densita', la sua pressione, la sua velocita';  • lo spazio viene rappresentato da una griglia di NxNxN elementi; • le variabili del fluido vengono calcolate sulle celle della griglia; • i valori calcolati sono i valori medi della variabili sulle celle; • questi valori sono calcolati risolvendo con una qualche approssimazione numerica le equazioni di conservazione: ...ad esempio, al primo ordine (scarsa accuratezza), in una dimensione, l'equazione di conservazione della massa si può approssimare alle differenze finite come:

14. A(i , j) Parallelizzazione di algoritmi fluidodinamici alle differenze finite E’ un problema “semplice”. Per loro intrinseca natura le equazioni fluidodinamiche sono locali. Pertanto la parallelizzazione del codice si può ottenere in modo efficiente attraverso una suddivisione a blocchi del dominio, tipo quella presentata per l’analisi delle immagini. Il codice risultante scala molto bene ed e’ molto ben bilanciato. Algoritmi di questo tipo si prestano molto bene alla parallelizzazione !!!

15. Esempio: Data la località della soluzione il processo N ha bisogno di conoscere SOLO la prima riga di dati dei processi N-1 e N+1 Si può dividere il dominio computazionale a “bande”. PE N+1 Ghost Region PE N Ghost Region PE N-1 L’array contenente la variabile viene definito con 2 righe in più, contenenti le “Ghost Regions”. La comunicazione coinvolge esclusivamente le ghost regions e i bordi. Il resto del calcolo è completamente locale.

16. L’algoritmo N-Body • Nel caso di simulazioni cosmologiche di sola materia oscura, si considera la sola massa, quindi le sole interazioni gravitazionali. In generale comunque la procedura puo' essere generalizzata a casi più articolati e complessi: • ogni particella ha massa M data; • la velocità di ogni particella viene calcolata in istanti di tempo successivi (discretizzazione del tempo) risolvendo le equazioni del moto con opportune tecniche alle differenze finite. • La forza gravitazionale che agisce su ciascuna particella viene calcolata attraverso un'opportuna approssimazione numerica.

17. Dinamica delle Particelle Approx. del primo ordine Equazioni del moto Approx. del secondo ordine: Leapfrog Approx. del secondo ordine: Lax Wendroff (two steps)

18. Approx. del primo ordine Calcolo delle Forze Metodo Particle Particle (PP): la forza su una particella è la risultante di tutte le interazioni con le altre particelle Vantaggi - accuratezza, risoluzione Limiti - numero di particelle (scala con N2) Metodi Tree:la forza è lo sviluppo di multipolo dovuto alla distribuzione delle particelle Vantaggi - risoluzione, scalabilità (N log(N)) Limiti - complessità della costruzione dell'albero

19. Metodo Particle Mesh la forza viene calcolata attraverso una doppia trasformata di Fourier sul campo di densità associato alla distribuzione di particelle ... Vantaggi - scalabilità, rapidità, semplicità dell'algoritmo Limiti - Risoluzione spaziale • discretizzo lo spazio con una griglia cubica si lato L, numero di punti per dimensione N e (quindi) risoluzione spaziale dx=L/N. • calcolo il campo di densità associato alla distribuzione di particelle distribuendo la massa di ciascuna particella sui punti griglia. • risolvo l'equazione di Poisson (che lega la densità al potenziale gravitazionale) attraverso una procedura FFT. In questo modo calcolo il campo gravitazionale sulla griglia. • calcolo la forza differenziando il potenziale e stimo la forza che agisce sulla particella con una procedura analoga alla 2).

20. Particle Mesh: distribuzione dei dati • Distibuzione dei dati tra i processori: In questo caso esistono due tipi di vettori • quelli che descrivono le particelle (posizione, velocita'); • quelli che descrivono i campi sulla griglia (potenziale gravitazionale, densita'). • questi ultimi vengono distribuiti a piani paralleli.Le particelle vengono suddivise tra i processori: NB: non c'e' correlazione tra la posizione delle particelle locali ad un processore, l'id del processore, e il piano della griglia locale a quel processore Questo è un problema essenziale per la parallelizzazione !!!

24. Particle Mesh: considerazioni. • in generale una particella sarà locale ad un processore diverso da quello che possiede i dati su griglia necessari a calcolarne la dinamica. • 2) il numero di particelle che "migrano" cambia continuamente e in modo disomogeneo man mano che la simulazione procede. • questo porta ad un algoritmo che in generale richiede molta comunicazione e forti sbilanciamenti nel carico di lavoro che ogni processore deve realizzare !!!

25. Particle Mesh: distribuzione del lavoro. La scelta delle modalità di distribuzione del lavoro è quindi particolarmente delicata per l'efficienza del codice. Nel caso più semplice, si può decidere se far svolgere il lavoro al processore che ha locale la particella o quello che ha locali i dati su griglia. Se lo svolge il PE che"possiede" la particelladevo fargli avere tutti i dati delle celle vicine: gli algoritmi standard di distribuzione della massa o di calcolo delle forze richiedono da8 a 27 celle. In questo modo comunque ogni processorecalcola la dinamica di NPART_TOT/NPES particelle (perfetto bilanciamento). Se il lavoro lo svolge il PE che"possiede" le cellequesto deve richiedere soltanto 6 informazioni (la posizione e la velocita') della particella. Tuttavia si potrebbe verificareche un solo PE esegue tutto il lavoromentre gli altri attendono inoperosi. La scelta quindi dipende dal bilancio tra queste due opposte esigenze: assicurare che la comunicazione sia minima e garantire un buon bilanciamento del carico di lavoro.

26. Particle Mesh: distribuzione del lavoro. Implementazione MPI • Nel esempio che segue, abbiamo invece implementato la seconda strada, che può portare a forti sbilanciamenti, ma ottimizza la comunicazione. L'implementazione è basata su MPI. Consta dei seguenti step: • Individuazione della posizione della particella sulla griglia • Creazione di liste di particelle associate ai vari processori • Comunicazione dei dati delle particelle tra i vari processori • Calcolo della densita' • Calcolo del potenziale • Calcolo della forza • Integrazione delle equazioni del moto • La strategia è quindi quella di sapere in anticipo quali particelle cadono in processori diversi da quelli a cui appartengono e comunicare le relative informazioni, in un unico processo di comunicazione, ai processori di competenza che a questo punto possono eseguire il calcolo in modo completamente locale

27. Particle Mesh: distribuzione del lavoro. Implementazione MPI Questo tipo di implementazione consente il controllo completo della comunicazione ma necessita di sostanziali modifiche rispetto alla versione sequenziale. Step 1: creo alcuni array per lo scambio di informazioni “generali” tra i processori allocate (count_pe_s(npes),STAT=error) allocate (count_pe_r(npes),STAT=error) count_pe_s=0 count_pe_r=0 ! ! find how many particles fall in each processor ! n_tot_s=0 do ipart=1,nparmax i3=int(x3(ipart)+0.5) if(i3.eq.0)i3=npes*nz ipe3=(i3-1)/nz if(ipe3.ne.mype)then count_pe_s(ipe3+1)=count_pe_s(ipe3+1)+1 n_tot_s=n_tot_s+1 endif enddo Ogni posizione in questi array contiene il numero di particelle che dovranno essere inviate/ricevute a/da un processore remoto. In questo loop l’array viene effettivamente riempito. Conto quante particelle in totale ogni processore invia

28. Particle Mesh: distribuzione del lavoro. Implementazione MPI Step 2: comunico il numero di particelle che i vari processori si scambiano. CALL MPI_Alltoall ( count_pe_s, 1, MPI_INTEGER, count_pe_r, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr ) …or, using send end recv… do i=1,npes if(i-1.ne.mype)then to_pe=i-1 from_pe=i-1 CALL MPI_ISend(count_pe_s(i),1,,to_pe,10,& MPI_COMM_WORLD,ierr) CALL MPI_IRecv( (i),1,,& from_pe,10,,status,ierr) endif enddo

29. Particle Mesh: distribuzione del lavoro. Implementazione MPI Step 3: Sapendo quante particelle devono venire comunicate, preparo i buffer con cui comunicare le rispettive posizioni e velocità. allocate (x1_aux_s(n_tot_s+1),STAT=error) allocate (x2_aux_s(n_tot_s+1),STAT=error) allocate (x3_aux_s(n_tot_s+1),STAT=error) allocate (v1_aux_s(n_tot_s+1),STAT=error) allocate (v2_aux_s(n_tot_s+1),STAT=error) allocate (v3_aux_s(n_tot_s+1),STAT=error) do ipart=1,nparmax i3=int(x3(ipart)+0.5) if(i3.eq.0)i3=npes*nz ipe3=(i3-1)/nz if(ipe3.ne.mype)then index0=point_pe_aux(ipe3+1) x1_aux_s(index0)=x1(ipart) x2_aux_s(index0)=x2(ipart) x3_aux_s(index0)=x3(ipart) v1_aux_s(index0)=v1(ipart) v2_aux_s(index0)=v2(ipart) v3_aux_s(index0)=v3(ipart) n_pepe(index0)=ipart point_pe_aux(ipe3+1)=point_pe_aux(ipe3+1)+1 endif enddo Step 4: Riempio i buffer

30. Particle Mesh: distribuzione del lavoro. Implementazione MPI Step 5: Allo stesso modo alloco i vettori che dovranno ricevere i dati. allocate (x1_aux_r(n_tot+1),STAT=error) allocate (x2_aux_r(n_tot+1),STAT=error) allocate (x3_aux_r(n_tot+1),STAT=error) allocate (v1_aux_r(n_tot+1),STAT=error) allocate (v2_aux_r(n_tot+1),STAT=error) allocate (v3_aux_r(n_tot+1),STAT=error) do i=1,npes if((i-1).ne.mype)then to_pe=i-1 from_pe=i-1 index0=point_pe(i) index1=point_pe_r(i) stride0=count_pe_s(i) stride1=count_pe_r(i) … CALL MPI_Irecv(x1_aux_r(index1),stride1,MPI_DOUBLE_PRECISION,& from_pe,20,MPI_COMM_WORLD,req(1),ierr) CALL MPI_Isend(x1_aux_s(index0),stride0,MPI_DOUBLE_PRECISION,& to_pe,20,MPI_COMM_WORLD,req(2),ierr) … CALL MPI_WAITALL(12,req,status_array,ierr) endif enddo

31. Particle Mesh: distribuzione del lavoro. Implementazione MPI Step 6: gli array di send non servono più. Per risparmiare spazio li elimino… deallocate (x1_aux_s) deallocate (x2_aux_s) deallocate (x3_aux_s) deallocate (v1_aux_s) deallocate (v2_aux_s) deallocate (v3_aux_s) Step 7: tutti i dati sono ora locali ai processori. Posso procedere con il calcolo come se fosse sequenziale. Calcolo la densità a partire dalle particelle call density(x1_aux_r,x2_aux_r,x3_aux_r,n_tot) call fields call gravity call nbody(n_tot,x1_aux_r,x2_aux_r,x3_aux_r,& v1_aux_r,v2_aux_r,v3_aux_r) Calcolo il campo gravitazionale Calcolo le forze gravitazionali Sposto le particelle

32. Particle Mesh: distribuzione del lavoro. Implementazione MPI Step 8: ripeto al contrario quanto fatto in precedenza, riportando le informazioni su posizioni e velocità delle particelle comunicate ai processori originari. Step 9: dealloco tutti gli array ausiliari necessari alla comunicazione. Questa operazione è utile a risparmiare memoria ma, soprattutto, è necessaria in quanto il numero di particelle da comunicare e quindi la dimensione degli array, cambiano ad ogni step evolutivo.

33. Particle Mesh: distribuzione del lavoro. Implementazione MPI 2 MPI 1 non consentendo accessi asincroni (one-side communication) alla memoria, necessita, in questo caso, di una “infrastruttura” piuttosto complessa. Con MPI 2 si può invece ottenere un’implementazione molto più snella, simile a quella Open MP - shared memory vista in precedenza. call MPI_FENCE(…) call MPI_WIN_CREATE(phi,8*N**3,MPI_DOUBLE_PRECISION,…) call MPI_FENCE(…) do j=1,npart ip=nint(x1(j)) trovo la posizione della particella sulla griglia jp=nint(x2(j)) kp=nint(x3(j)) If(la particella è “locale”)then forza_media=...    else call MPI_LOCK(MPI_LOCK_SHARED,…) call MPI_GET(dati necessari per il calcolo della forza) call MPI_UNLOCK(…) forza_media=… endif v_new = v_old + forza_media * dtaggiorno la velocita' x_new = x+old + v_new *dt aggiorno la posizione enddo call MPI_WIN_FREE(…)

34. Tree Code In questo caso NON viene introdotta alcuna griglia computazionale... La forza sulla particella i-esima viene approssimata come somma dei contributi delle particelle "vicine": ... e di quell lontane, attraverso un'approssimazione di multipolo del potenziale gravitazionale: ove con il pedice cm viene indicato il centro di massa dell'intero set di particelle e con M la massa totale. Tipicamente si blocca lo sviluppo in serie all'ordine di quadrupolo

35. Tree Code Quali i vantaggi ??? Il multipolo viene calcolato non sulle particelle, ma su un sottoinsieme di celle. Le celle vengono calcolate secondo una procedura ad albero: step 1:definisco un angolo critico A step 2:divido il box computazionale in 8 celle ciascuna di dimensione l. Calcolo il centro di massa di ciascuna cella e calcolo la distanza d tra la particella e il centro di massa. Se l/d > Ala cella è lontana e di essa si calcola il quadrupolo, se l/d < Asuddivido ulteriormente la cella step 3:ripeto la procedura fino a quando rimango con cellette vuote o contenenti un'unica particella. A queste applico la legge di Newton.

36. Tree Code: implementazione parallela Viene attribuito lo stesso numero di particelle ad ogni processore. Utilizzando l'albero le particelle vengono redistribuite dinamicamente tra i processori durante l'evoluzione in modo che particelle fisicamente vicine siano anche nello stesso processore) • Due tipi di vettori fondamentali: • Particelle (posizione, velocità, e massa) • Tree (posizione, massa del baricentro) • La distribuzione dei dati è fatta in modo da garantire una distribuzione omogenea della memoria. Il lavoro è distribuito equamente tra le particelle (load balancing) Le celle dell'albero vengono numerate progressivamente dalla più grande alle più piccole. I livelli più grossolani vengono distribuiti in modo ciclico, in modo da evitare il problema di accesso a risorse critiche: se i livelli più grossolani fossero tutti sullo stesso processore tutti andrebbero ad accedere alla memoria di quel processore. I livelli più fini sono invece locali al processore che possiede il "padre".

37. Tree Code: implementazione parallela Problema fondamentale del tree code. Calcolo della forza: ogni particella deve recuperare i dati remoti di tutte le particelle "vicine" che non stanno nello stesso processore tutti i quadrupoli delle celle "lontane" Calcolo del tree: per calcolare la componenente di quadrupolo della cella può essere necessario accedere a particelle "remote" (nel senso della memoria) La comunicazione è estremamente pesante vengono utilizzate one side comunication (latenza minima) basate su: shmem    T3E Origin 3000 LAPISP3 Implementazione MPI 2 è in corso

38. Tree Code: ottimizzazione • 1 Grouping • Particelle molto vicine possono risentire approssimativamente della stessa interazione di quadrupolo da celle lontane.Pertanto è inutile ricalcolarla per ogniuna di queste particelle. • E' possibile quindi: • definire due distanze caratteristiche R e D; • partendo da posizioni opportune (zone particolarmente dense) raggruppare le particelle che stanno in una sfera di raggio R; • calcolare la forza su queste particelle come: Ftot = Fnear + Ffar,   ove:Fnearviene calcolata in modo "standard" per particelle e celle che distano meno di D,Ffarviene calcolata una volta per tutte per particelle che distano più di D.

39. Tree Code: ottimizzazione 2 Dynamic Load Balancing A causa dei diversi tempi di calcolo delle forze per ciascuna particella, un processore può terminare prima degli altri. In questo caso (grazie alle one-side comunications, che non presumono comunicazione) può "aiutare" un'altro processore a completare il suo carico di lavoro. 3 Data Buffering Spesso gli stessi dati, relativi a particelle o celle remote, servono per il calcolo delle forza che agisce su particelle diverse. Pertanto è stato implementato un meccanismo di "cache". La memoria del processore non allocata, viene riempita (per quanto possibile) con dati prelevati in precedenza su altri processori, che a quel punto divengono locali.