Dipartimento di Fisica e Astronomia - Indiceinsieme di equazioni differenziali ordinarie e le forze...

Indice 1. Introduzione ...............................................................................................4

2. Gravitazione e forze idrodinamiche ..........................................................8

2.1 Introduzione .........................................................................................8

2.2 Smoothed Particle Hydrodynamics .....................................................9

2.2.1 Introduzione ..................................................................................9

2.2.2 Proprietà .....................................................................................10

2.2.3 Calcolo della densità...................................................................12

2.2.4 Calcolo della pressione ..............................................................13

2.2.5 Calcolo delle forze ......................................................................14

2.2.6 Lunghezza di smussamento ........................................................15

2.3 Gravitazione ......................................................................................16

2.3.1 Introduzione ................................................................................16

2.3.2 Algoritmo classico ......................................................................17

2.3.3 Algoritmo di Barnes-Hut............................................................17

2.3.3.1 Introduzione .........................................................................17

2.3.3.2 Costruzione dell’oct-tree .....................................................18

2.3.3.3 Costruzione dei nodi ed inserimento delle particelle

nell’albero ........................................................................................19

2.3.3.4 Calcolo del centro di massa per ogni nodo..........................20

3. Reimpostazione del codice ......................................................................21

INDICE 2

3.1 Introduzione .......................................................................................21

3.2 Descrizione logica .............................................................................22

3.3 Modifiche strutturali ..........................................................................24

3.4 Modifiche dell’interfaccia utente ......................................................25

3.4.1 Miglioramento della linea di comando .......................................26

3.4.2 Miglioramento dell’output a video .............................................27

3.4.3 Miglioramento del file di configurazione ...................................28

3.4.4 Aggiunta di file di log e statistici ...............................................29

3.5 Modifiche funzionali .........................................................................30

3.5.1 Arresto forzato di una simulazione .............................................30

3.5.2 Riavvio di una simulazione ........................................................31

3.5.3 Supporto di file little e big endian..............................................32

3.5.4 Scelta manuale del tempo di creazione di ogni snapshot ...........33

3.6 Strutture dati e variabili principali.....................................................34

3.7 Implementazione dell’algoritmo di Barnes-Hut ................................37

3.7.1 Introduzione ................................................................................37

3.7.2 Implementazione ........................................................................38

3.7.3 Aggiornamenti dinamici dell’albero...........................................40

3.7.3.1 Complessità .........................................................................41

3.7.4 Riallocazione dinamica dei nodi ................................................41

3.8 Nuovo criterio per l’apertura dell’oct-tree ........................................42

3.9 Supporto delle coordinate comoventi................................................45

3.10 Supporto delle condizioni periodiche al contorno ...........................46

3.11 Integrazione temporale ....................................................................47

3.11.1 Integratore temporale per una singola particella ......................47

3.11.2 Integratore temporale per sistemi N-Body...............................49

INDICE 3

3.12 Ricerca dei vicini .............................................................................51

3.13 Moduli utilizzati ..............................................................................54

3.14 Albero delle chiamate ......................................................................64

3.15 Valutazione dell’efficienza ..............................................................71

4. L’utility StoB...........................................................................................77

4.1 Introduzione .......................................................................................77

4.2 Caratteristiche principali ...................................................................78

4.3 Uso di StoB dalla linea di comando ..................................................80

4.4 Generazione di immagini in formato bmp.........................................81

4.5 Determinazione automatica o manuale del dominio .........................82

4.6 Rappresentazione delle particelle con colori e simboli diversi .........83

4.7 Supporto delle rotazioni e della prospettiva ......................................84

4.8 Immagini “quasi” fotografiche ..........................................................85

4.9 Il linguaggio di script.........................................................................86

5. Conclusioni..............................................................................................90

A. Glossario .................................................................................................93

B. Il file dei parametri .................................................................................99

C. Alberi rosso-neri ...................................................................................113

D. Il modulo binary_m ..............................................................................117

Bibliografia ................................................................................................123

Indice delle figure ......................................................................................126

Indice analitico ..........................................................................................127

CAPITOLO 1. INTRODUZIONE 4

1. Introduzione Le simulazioni numeriche tridimensionali sono diventate indispensabili nella cosmologia, infatti numerosi problemi astrofisici e cosmologici richiedono la simulazione dell’evoluzione temporale di un sistema fisico. Si pensi, ad esempio, allo studio della formazione dei aloni di materia oscura, degli ammassi di galassie e l’interazione tra galassie isolate, etc. Senza le tecniche numeriche, gli immensi progressi fatti in questo campo non sarebbero stati possibili. La complessità di questi problemi, è tale che per contenere i tempi di calcolo (anche su elaboratori molto veloci) è necessario adottare delle approssimazioni di tipo numerico. La creazione di algoritmi sempre più veloci rimane l’interesse principale per sfruttare appieno i vantaggi offerti dai moderni elaboratori. In questo tipo di simulazioni sono molto importanti gli algoritmi per il calcolo della forza gravitazionale tra le particelle. Usando l’approccio diretto, ovvero calcolando la forza gravitazionale agente su di una particella come la somma delle forze che agiscono su di essa dovute alle particelle rimanenti, la complessità del metodo è molto elevata. Infatti il costo computazionale, seguendo questa implementazione, risulta essere pari a )( 2NO ed è quindi impensabile usare questa tecnica in simulazioni con un gran numero di corpi, anche su computer molto veloci. Nel tempo sono state sviluppate numerose tecniche per calcolare campi gravitazionali su larga scala.


Noi abbiamo usato un algoritmo ad albero ideato da Barnes e Hut nel 1986, che permette di abbassare la complessità del calcolo delle forze gravitazionali a )log( NNO . Nelle nostre simulazioni abbiamo fondamentalmente due tipi di particelle: le particelle barioniche e le particelle di materia oscura. Mentre le seconde sono soggette solamente a forze di tipo gravitazionale, le prime sono soggette anche a campi di tipo idrodinamico, come la temperatura e la pressione. Per calcolare questi campi, abbiamo impiegato la tecnica numerica chiamata SPH (Smoothed Particles Hydrodynamics), supponendo che tutta la materia barionica sia gas. Questo metodo ben si adatta ad ogni tipo di geometria e la sua natura Lagrangiana permette un cambio locale di risoluzione che automaticamente segue la locale densità di massa. Questa caratteristica è molto importante, poiché permette di risparmiare tempo di elaborazione concentrando i calcoli laddove si hanno alte concentrazioni di gas.

Figura 1: condizioni iniziali della simulazione

Campo di densità delle particelle al momento iniziale; a sinistra la materia oscura, a destra il gas.


In questa tesi illustreremo tutte queste tecniche, descrivendo le modifiche che sono state fatte al codice da cui siamo partiti [D. Dal Borgo, 2000]. Dapprima verranno esposti i concetti fisici legati al problema, seguiti da una breve ma completa introduzione all’algoritmo di Barnes-Hut. Nel capitolo 3 verrà presentato il programma che abbiamo chiamato “inSPHector”: verranno descritti tutti gli accorgimenti che abbiamo implementato per aumentarne l’efficienza, la robustezza e la leggibilità. Una serie di ottimizzazioni ci hanno permesso di migliorare notevolmente le prestazioni del programma:

• abbiamo implementato un BH-tree per ogni tipo di particella; • il calcolo dei vicini nel metodo SPH ora è effettuato usando un

algoritmo che percorre il BH-tree stesso; • è stato ottimizzato l’integratore temporale; • nuovo criterio di apertura dell’oct-tree.

È stato implementato un supporto delle coordinate comoventi e delle condizioni periodiche al contorno. Abbiamo reso il codice:

• modulare; • leggibile e ben strutturato: sono state ampiamente usate le

“strutture” del Fortran 90 (si veda il glossario); • adattivo: la memoria viene allocata dinamicamente.

Abbiamo rivolto particolare attenzione all’interazione dell’utente con l’applicativo; abbiamo riorganizzato il codice originale in modo che il programma:

• sia semplice da usare; • sia protetto da eventuali errori commessi dall’utente; • supporti files sia in formato big endian che in formato little endian; • preveda numerose opzioni.


Figura 2: risultato di una simulazione idrodinamica

Distribuzione finale dei campi di densità della materia oscura (in alto a sinistra) e del gas (in alto a destra).

Infine, nel capitolo 4, verrà presentata l’utility StoB (Snapshot to bitmap) indispensabile per il debugging del programma.

CAPITOLO 2. GRAVITAZIONE E FORZE IDRODINAMICHE 8

2. Gravitazione e forze idrodinamiche

2.1 Introduzione In cosmologia, problemi come la formazione e l’evoluzione di una galassia o di un ammasso di galassie richiedono sia la risoluzione di equazioni che regolano l’evoluzione di un sistema soggetto a forze gravitazionali sia il calcolo delle forze idrodinamiche. Nelle nostre simulazioni vi sono due tipi di particelle: le particelle barioniche e le particelle di materia oscura. Le prime rappresentano la materia ordinaria, ovvero quella visibile, come l’idrogeno, l’elio, il ferro, etc. Le seconde, invece, rappresentano quella parte della materia a noi non visibile la cui esistenza è stata provata solo per via indiretta. Questo tipo di materia potrebbe essere costituita da particelle fisiche che hanno la caratteristica di interagire solo debolmente col resto: in pratica sono soggette solamente a forze gravitazionali. I dati osservativi suggeriscono che ben il novanta percento dell’universo sia costituito da questo tipo di materia. Il restante dieci percento è costituito da materia barionica, essenzialente idrogeno, con tracce di altri elementi come elio, carbonio, ferro, etc.


In un sistema da evolvere, le particelle sono caratterizzate da una serie di valori alcuni dei quali sono dipendenti dal tempo come l’accelerazione e la velocità. Le particelle con accelerazioni a più alte devono essere mosse più di frequente rispetto a quelle con accelerazione minore, in modo da ottenere una maggiore accuratezza dei risultati. Si può scegliere

a1

∝∆t oppure a

1∝∆t (si veda paragrafo 3.11 per maggiori dettagli).

2.2 Smoothed Particle Hydrodynamics

2.2.1 Introduzione Per evolvere correttamente un sistema di particelle, si deve prevedere nel codice un supporto per il calcolo delle forze di tipo idrodinamico. Per fare questo abbiamo usato il metodo chiamato Smoothed Particle Hydrodynamics (che abbrevieremo con SPH). Esso è una potente tecnica Lagrangiana nata per la risoluzione di equazioni dell’idrodinamica e rappresenta una valida alternativa ai metodi euleriani che necessitano di una griglia. Poiché nello spazio il gas (idrogeno per la precisione) costituisce quasi la totalità della materia visibile, d’ora in avanti useremo interscambiabilmente i termini “gas” e “materia barionica”. La tecnica SPH rappresenta il gas come un insieme di particelle ed il suo sistema di equazioni regola l’evoluzione ed il comportamento di ciascuna di esse. La sua natura lagrangiana garantisce una maggior risoluzione spaziale


nelle zone più dense, ovvero dove ci sono più particelle, mentre nelle zone meno dense la stima è meno precisa. La natura particellare del metodo SPH semplifica notevolmente il calcolo numerico delle derivate spaziali, che si riducono a interpolazioni sulle posizioni delle particelle. Le equazioni dell’idrodinamica diventano un insieme di equazioni differenziali ordinarie e le forze tra particelle di gas sono rappresentate come interazioni tra coppie di particelle.

2.2.2 Proprietà La tecnica SPH, per descrivere un fluido, si basa su due assunzioni fondamentali: 1. le proprietà f del fluido sono stimate in ogni punto prendendo una media

pesata in un intorno del punto stesso:

'),'()'()( rrrrr dhWfffaa

−== ∫ (2.1)

dove h è una lunghezza caratteristica, che rappresenta, in un certo senso, le dimensioni dell’elemento di fluido, r è il vettore posizione della particella e W è una funzione di peso detta “kernel”. W deve avere le seguenti proprietà: • 1'),'( =−∫ rrr dhW (2.2)

• )'(),'(lim0

rrrr −=−→

δhWh

(2.3)

dove )(xδ è la funzione di Dirac. Diciamo che f è la quantità

“smussata” (smoothed).


2. Un fluido continuo è schematizzabile come un numero finito N di particelle, ognuna delle quali è libera di muoversi sotto l’influenza delle altre; si può approssimare l’integrale (2.1) con:

( ) ∑=

−=N

jjij

j

ji hWf

mf

1

),()( rrrrρ

, (2.4)

dove jm è la massa della j-esima particella e jρ la densità nella sua

posizione. Se però imponiamo 0),( =− hW ji rr per mhji >− rr , in cui m è

un intero opportunamente scelto, quando si considera la particella

ip , solo le particelle Ip j ∈ , con

mhpppI jij <−= | ,

contribuiscono alle sommatorie. Una scelta opportuna di m consente di avere ff

a≅ e contemporaneamente di ridurre il

costo per calcolare f .

Sviluppando (2.1) in serie di Taylor nel punto rr =' si ottiene:

).()()()( 322 hOhfcff i +∇+= rr (2.5)

Approssimando f con f si ottiene un errore )( 2hO .


2.2.3 Calcolo della densità Col metodo SPH, per calcolare la densità di una particella i, usiamo:

∑=

=N

jiijji hWm

1

),(rρ (2.6)

dove jiij rrr −≡ e im è la sua massa.

La funzione kernel che abbiamo preso in considerazione è una spline cubica ed è stata introdotta da Monaghan e Lattazio nel 1995:

≤<

−

≤≤

+

−

=

altrimentihr

hr

hr

hr

hr

ChrW

0

121

12

21

0661

),(3

32

(2.7)

dove la costante C dipende da h:

3

1h

Cπ

= . (2.8)


Figura 3: funzione kernel usata nel codice

2.2.4 Calcolo della pressione Per calcolare la pressione della particella i usiamo:

iii uP ργ )1( −= (2.9)

in cui γ è un coefficiente chiamato “esponente adiabatico del gas” (nel

nostro caso 35

=γ ), mentre iu rappresenta la sua energia interna e può

essere ricavato dalla temperatura della particella in esame:


p

ii m

kTu

µγ )1( −= (2.10)

dove k è la costante di Boltzmann, µ è il peso molecolare medio e pm è la

massa del protone.

2.2.5 Calcolo delle forze Sia v il vettore velocità. L’equazione di Eulero:

ρP

t∇

−=∇⋅+∂∂

vvv

)(

lega la pressione P alla densità ρ . Si dimostra (si veda [Dal Borgo, 2000]) che l’accelerazione della particella i di gas si può calcolare con:

∇+∇

∏++−=+

∇−= ∑ ),(

21

),(21~

22 jijiiijij

ijj

j

i

ij

visci

i

gasi hWhW

PPm

Prraa

ρρρ

. (2.11)

La sua variazione di energia interna è:

∇+∇−

∏++= ∑ ),(

21

),(21

)(~

21

22 jijiiijij

jiijj

j

i

ij

i hWhWPP

mdtdu

rrvvρρ

.

(2.12)


Il termine ij∏~

rappresenta la viscosità artificiale e vale:

( ) ijjiij ff ∏+=∏21~

(2.13)

con

<⋅

+−=∏

altrimenti

sec

ijijij

ijijij

ij

0

02 2

rvρ

αµµα, (2.14)

e ( )

( ) ( )ii

iif

vv

v

×∇+⋅∇

⋅∇= (2.15)

ed infine

( )( )22

ijji

jijiijij

h

h

εµ

+−

−−=

rr

rrvv (2.16)

Nelle equazioni (2.11) e (2.12) una particella i interagisce con una particella j solo se iij h<r o jij h<r , dove h viene chiamata lunghezza di

smussamento.

2.2.6 Lunghezza di smussamento Essa rappresenta la dimensione del più piccolo elemento di fluido che SPH può risolvere. Da un punto di vista pratico, la lunghezza h è legata al numero di particelle che contribuiscono alla stima nelle equazioni idrodinamiche e, nel nostro


codice, è variabile in modo che ogni particella interagisca sempre con altre n<N particelle (n impostato nel file dei parametri alla voce SPHNeighbours). È stato dimostrato, infatti, che se il numero vicini è costante si ottiene una migliore accuratezza dei risultati e una migliore conservazione dell’energia.

2.3 Gravitazione

2.3.1 Introduzione L’equazione nota come legge di Newton, che regola la forza gravitazionale tra un corpo con massa 1m e un corpo con massa 2m , distanti tra loro 12r è la seguente:

312

21121

rr

mmGF =

r, (2.17)

G è la costante di gravitazione universale, pari a 6.672e-8 123 sec −− gcm . Dati N corpi, la forza che N-1 corpi esercitano sul rimanente è:

3

ij

ji

jiijj

mmG

rrF ∑

≠

= i=1,...,N (2.18)

Questa equazione afferma che la forza agente su una particella, dipende dalla posizione e dalle masse di tutte le rimanenti.


Essa è additiva: la forza che agisce su un corpo è pari alla somma di tutte le singole forze esercitate dai rimanenti corpi del sistema.

2.3.2 Algoritmo classico È facile scrivere un algoritmo che, date N particelle, calcola la forza gravitazionale agente su ognuna di esse. Poiché per calcolare la forza agente su una particella ci servono N-1 termini, in totale ce ne sono )1( −NN . Da una nota proprietà della legge di gravitazione abbiamo che:

jiij FF −=

Quindi, in totale, il calcolo della forza costa 2

)1( −NN operazioni.

L’algoritmo ha quindi complessità )( 2NO . Sono stati inventati numerosi algoritmi e strutture dati che cercano di abbassare questa complessità, che è proibitiva nelle simulazioni con un gran numero di particelle. Uno di questi algoritmi è quello di Barnes-Hut, che descriveremo ora.

2.3.3 Algoritmo di Barnes-Hut

2.3.3.1 Introduzione

Questo algoritmo, ideato nel 1986, è un algoritmo di tipo gerarchico e adattivo che ben si presta al calcolo delle forze gravitazionali. L’idea di base è quella di rendere più veloce possibile il calcolo delle forze gravitazionali, usando opportune approssimazioni.


Si approssima la formula (2.18) con la seguente:

3

,

,3

cmi

icmi

jiij

jiiji

MmG

mmG

rr

rrF ≈= ∑

≠

i=1,.....N (2.19)

Essa asserisce che per ogni singolo corpo i, il contributo di forza esercitato su di esso da un agglomerato composto da knn ........1 elementi

sufficientemente lontani, può essere approssimato con un singolo punto p posto nel centro di massa di knn ........1 .

Risulta naturale, quindi, pensare di calcolare la forza esercitata su una particella, come la somma vettoriale della forza calcolata con la formula (2.19), per le particelle più lontane di una certo valore m, e la forza calcolata usando la formula (2.17), per le altre. Questo metodo abbassa notevolmente i tempi di calcolo e raggiunge maggiore o minore precisione in base ad m. In pratica più m è piccolo, meno accurati sono i risultati e più veloce è il calcolo. L’equazione (2.19) non è valida se la particella su cui si intende calcolare la forza è dentro una sfera che contiene tutti i corpi che si approssimano con il centro di massa. L’algoritmo di Barnes-Hut si basa su una suddivisione ricorsiva del dominio della simulazione, basandosi sulla densità relativa della regione considerata. A tale scopo è stata introdotta una struttura dati ad albero, chiamata oct-tree.

2.3.3.2 Costruzione dell’oct-tree La costruzione dell’oct-tree avviene in due fasi:

• costruzione dei nodi ed inserimento di ogni particella nell’albero; • calcolo del centro di massa per ogni nodo.


2.3.3.3 Costruzione dei nodi ed inserimento delle particelle nell’albero

Inizialmente si racchiude il sistema in una cella madre cubica contenente tutte le particelle del sistema. Le particelle vengono inserite una alla volta e in caso di appartenenza allo stesso sottocubo, quest’ultimo viene suddiviso in otto parti. Questo procedimento crea quindi, ad ogni passo, una suddivisione dello spazio, creando così un albero gerarchico, dove ogni nodo deve “contenere” almeno una particella. L’albero è costituito da particellenodi nn + elementi.

Ogni nodo può avere al massimo otto possibili discendenti. Per distinguere i nodi dalle particelle, abbiamo enumerato da 1 a particellen le particelle e da

1+particellen a particellenodi nn + i nodi.

Questo è particolarmente importante quando si percorre l’albero, per decidere se il nodo che si sta esaminando è una cella ancora suddivisibile. I nodi che contengono delle sottocelle non vuote hanno come massa la somma delle masse delle sottocelle e come centro di massa il centro di massa delle particelle contenute nelle sottocelle.

Figura 4: esempio di oct-tree


La complessità computazionale della costruzione di un oct-tree è direttamente proporzionale alla densità delle particelle interne allo spazio di simulazione. Se le particelle sono uniformemente distribuite, si ottiene un albero bilanciato in cui le foglie (particelle) sono posizionate tutte allo stesso livello. La complessità risultante è quindi )log( 8 NNO .

Nel caso di distribuzione non uniforme la complessità può essere considerata pseudo-lineare )(NbO , con bN ≤2log , poiché esistono al più

b2 particelle distinte, dove b è il numero di bit usati per rappresentare la posizione delle particelle.

2.3.3.4 Calcolo del centro di massa per ogni nodo Costruito l’oct-tree, si associa a ciascun nodo il centro di massa del sistema particelle-nodo in esso contenuto mediante visita bottom-up dell’albero. Se supponiamo che l’albero sia bilanciato, la complessità computazionale è

)log( 8 NNO , poiché è necessario compiere la visita di un albero di N

foglie.

CAPITOLO 3. REIMPOSTAZIONE DEL CODICE 21

3. Reimpostazione del codice

3.1 Introduzione Il codice da noi rielaborato è l’evoluzione di un codice scritto da più persone in tempi diversi. La versione originale del programma è stata scritta negli anni ’80 da W.Benz. In essa venivano risolte le equazioni idrodinamiche per particelle immerse in un potenziale gravitazionale esterno in un contesto non cosmologico. Come integratore temporale veniva usato il metodo Runge-Kutta: i timestep individuali non erano arbitrari per ciascuna particella e dovevano assumere valori multipli di due. Nel 1991, Navarro e Benz modificarono il codice introducendo l’autogravità e la possibilità di trattare particelle non collisionali, ovvero di tutte quelle particelle che non seguono le leggi idrodinamiche, come ad esempio la materia oscura. Essi, inoltre, aggiunsero anche la trattazione di alcuni processi fisici importanti in astrofisica, come il raffreddamento del gas, la formazione stellare e gli effetti energetici dovuti alle supernovae. Nel 1997, Wendy Groom sostituì lo schema di integrazione Runge-Kutta con un “central leap-frog” e generalizzò la scelta del timestep permettendo di associare alle particelle un intervallo di tempo arbitrario. Infine, nel 2000, Davide Dal Borgo ha trasformato buona parte del codice, scritto originariamente in Fortran 77, nel più moderno Fortran 90.


Ha inoltre sostituito l’algoritmo di Appeal per il calcolo delle forze gravitazionali, con l’algoritmo di Barnes-Hut (1986), rendendo così il codice circa tre volte più veloce. Questo codice però soffre di almeno due fondamentali limitazioni: non supporta né le coordinate comoventi, né le condizioni periodiche al contorno. Con questa tesi abbiamo rimosso queste ed altre limitazioni, ottenendo così un codice moderno ed efficiente. Le modifiche al codice sono di tre tipi: modifiche strutturali, funzionali e dell’interfaccia utente. Esse verranno discusse in dettaglio nei prossimi paragrafi.

3.2 Descrizione logica Una simulazione è caratterizzata da un tempo iniziale e uno finale. Come già accennato al paragrafo 2.1, le particelle con accelerazioni più alte devono essere mosse più di frequente rispetto a quelle con accelerazione minore in modo da ottenere una migliore accuratezza dei risultati. Ad ogni particella viene associata un’informazione che indica il suo tempo corrente e un’informazione t∆ che indica dopo quanto tempo, rispetto al tempo corrente, essa deve essere mossa. Una volta letti tutti i valori necessari per una corretta simulazione è necessario calcolare le accelerazioni delle particelle. Per fare questo si devono calcolare le forze idrodinamiche e quelle gravitazionali. Nel passo successivo si calcolano i t∆ per tutte le particelle e si spostano solamente quelle che soddisfano il criterio descritto al paragrafo 3.11.2. Si ricalcolano le accelerazioni per le particelle spostate ed il ciclo si ripete.


Il programma termina quando il tempo della simulazione è maggiore o uguale al tempo massimo impostato nel file dei parametri. Ci sono altri modi per terminare una simulazione; essi sono descritti al paragrafo 3.5.1.

Calcolo accelerazione particelle

Gravitazione

File dei parametri

Condizioni iniziali

E' stato trovato il file "stop"?

Ricerca delle particelle da muovere percorrendo

l'albero rosso-nero

Calcolo accelerazione

particelle

Avanzamento particelle

Calcolo dei nuovi timestep per le

particelle mosse

E' stato superatoil tempo tra due

file di restart?

E' finito il tempo CPU a disposizione?

E' stato superato il tempo finale della simulazione?

No

No

Start

file di riavvio

file di riavvio

E' stato superatoil tempo tra due

snapshot?snapshot

Yes

Yes

Yes

No

Stop

Yes

file di riavvio snapshot

Yes

Aggiornamento oct-tree

Visita (top-down) dell'oct-tree.Termina la visita di un ramo se viene soddisfatto il criterio di apertura. Calcolo delle forze

gravitazionali.

SPH

Ricerca dei vicini e visita dell'oct-tree

Calcolo forze idrodinamiche

Calcolo forze totali

No

Devono essere scritte le statistiche per l'energia?

statistiche

YesCalcolo potenziale

No

No

Figura 5: logica del programma


Durante una simulazione la nostra applicazione scrive dei file che chiameremo snapshot, i quali rappresentano delle vere e proprie fotografie del sistema simulato e permettono, inoltre, il riavvio di una simulazione.

3.3 Modifiche strutturali Le modifiche strutturali sono tutti quei cambiamenti al codice atti ad aumentarne la leggibilità, la manutenibilità, l’omogeneità e la modularità. Le principali modifiche di questo tipo effettuate nel codice sono:

• il codice è scritto in standard ANSI-Fortran 90, sfruttando a fondo tutte le caratteristiche che questo linguaggio ci offre. Uniche eccezioni a questa regola sono costituite dall’uso di:

- sizeof(var): questa funzione risulta molto comoda, poiché restituisce la grandezza in byte della variabile che le si passa come parametro. Per rendere il codice standard, ma meno elegante e sicuro, si possono sostituire tutte le ricorrenze di questa funzione con un valore numerico pari al numero in byte della variabile var.

- \: questo carattere, usato assieme ad una write o ad una print, permette di stampare a video una stringa senza andare a capo riga. Es: WRITE (*,'(\,A)'),"..."

Entrambe le funzionalità sono comunemente implementate in tutti i moderni compilatori Fortran 90 e quindi, a nostro avviso, l’uso di esse non è limitativo, visto che possono essere facilmente sostituite (a scapito della leggibilità del codice, della sua manutenzione e della stampa a video).

• uso dei moduli. In questa nuova versione del programma sono stati usati massicciamente i moduli, visto che essi ci permettono di rendere il codice modulare e protetto. Il Fortran 90, tramite il concetto di modulo (si veda il glossario), permette di definire l’ambito d’uso di una variabile, funzione o


subroutine, attraverso l’uso di due parole chiave: public e private. Se una variabile, funzione o subroutine è dichiarata private, allora è visibile solo all’interno del modulo in cui è stata definita. Se invece è dichiarata public (o se non viene specificato alcun attributo), allora è visibile anche al di fuori del modulo. Questo meccanismo ci permette di vedere un modulo come una specie di “scatola chiusa”, in cui un utente può usare solo le variabili, le funzioni e le subroutine previste da chi ha scritto il codice, senza correre il rischio di “side-effects”.

• Sono stati usati i type (record) del Fortran 90 (si veda il glossario), al posto dell’uso degli array, laddove il codice si prestava. Allo stato attuale, ad esempio, ogni particella e ogni nodo del oct-tree, sono rappresentati da un record.

• Uso massiccio della memoria dinamica. Viene infatti allocata a run-time solo la memoria strettamente necessaria all’esecuzione del programma. Unica eccezione è rappresentata dalla modalità con cui è stato implementato l’oct-tree.

3.4 Modifiche dell’interfaccia utente Sono sia le modifiche al codice atte a rendere più leggibili e comprensibili gli output del programma, sia le modifiche atte a rendere più semplice l’uso del programma stesso. Le aggiunte principali sono:

• miglioramento della facilità d’uso del programma da linea di comando;

• miglioramento dell’output a video; • miglioramento del file dei parametri di configurazione; • aggiunta di numerosi file e informazioni statistiche.


Essi verranno descritti nel dettaglio qui di seguito.

3.4.1 Miglioramento della linea di comando È stato migliorata l’invocazione del programma dalla linea di comando, con l’aggiunta di numerose opzioni e di un semplice help. Eseguendo il programma, senza nessun parametro, viene visualizzato il seguente output a video: Usage: inSPHector [-i] [[-g|-r[:<Snapshot>]] <ParameterFile>] [--info] [[--generate|--restart[:<Snapshot>]] <ParameterFile>] -i (--info) Version and current precompiled options -g (--generate) <ParameterFile> Generate a template parameter file <ParameterFile> -r(--restart) <ParameterFile> Restart the simulation from restart file -r:(--restart:)<NUM> <ParameterFile> Restart the simulation from a snapshot with suffix <NUM> and prefix specified in <Parameterfile> (item SnapshotPrefix). <NUM> must be a

number. If the snapshot time exceeds MaximumTime, MaximumTime is reset to the snapshot time.

È possibile effettuare operazioni come la partenza ed il riavvio di una simulazione precedentemente interrotta. È possibile generare un file di configurazione test, e, infine, è possibile visualizzare a video alcune informazioni sul programma come versione, data di creazione, etc.


3.4.2 Miglioramento dell’output a video Partendo dal presupposto che l’utente deve sapere in ogni istante che operazioni sta compiendo il programma in esecuzione, abbiamo notevolmente migliorato l’output a video. Le informazioni di tipo quantitativo adesso vengono scritte su file. Si è prestata, inoltre, particolare attenzione alla formattazione dell’output a video ottenendo come risultato una visualizzazione pulita e ordinata. Un esempio di output è il seguente: . . . Start Ngb-tree construction...done Ngb-tree. Start force computation...done. Start potential computation...done. Writing snapshot file 'snapshot_000' at time 0.0000000...done. Begin Timestep 0, Time: 0.0251789, Timestep: 0.0251789 Start force computation...done. Start potential computation...done. . . .


3.4.3 Miglioramento del file di configurazione Il file di configurazione (o file dei parametri) è un file di testo che viene letto dal programma all’inizio di ogni simulazione. In esso devono essere specificati tutti i parametri cosmologici necessari per definire il sistema fisico che si vuole far evolvere e tutti i parametri numerici che determinano la dinamica del sistema. In questo file, inoltre, vengono specificati i nomi dei file di output, la directory di destinazione dei file di output, etc. Rispetto alla versione precedente, sono stati aggiunti diversi parametri per rendere il programma più flessibile, stando attenti però a non complicare troppo il file di configurazione. È possibile chiedere al programma di generare un file di configurazione test, di nome “nome file”, eseguendo il comando:

inSPHector.exe -g <nome file>

Questo file serve solamente come traccia. L’utente, una volta generato, deve modificarlo per renderlo congruo alle proprie esigenze. Ogni parametro in questo file è del tipo: nomeparametro <valore> Non è necessario che i parametri siano in un particolare ordine. Devono essere, però, tutti presenti. Il programma effettua un syntax checking del file prima di far partire la simulazione e avverte l’utente se riscontra degli errori e/o delle incongruenze. Non è possibile, inserire parametri non previsti, come non è possibile inserire un parametro più di una volta.


Il carattere di commento è il carattere “%”. Tutto ciò che lo segue, fino alla fine della riga, non viene preso in considerazione. Sono previste, inoltre, due istruzioni:

• %print <string> stampa a video <string> durante il parsing del file; • %sleep <sec> aspetta “sec” secondi durante il parsing del file.

3.4.4 Aggiunta di file di log e statistici All’inizio di ogni simulazione vengono creati tre file di testo che conterranno:

• tutti i passi seguiti dal programma (simile all’output a video); • statistiche sull’energia del sistema; • profiling e statistiche sui tempi di esecuzione delle principali parti

del programma. Se durante la compilazione di inSPHector viene utilizzata l’opzione di precompilazione DIAG, all’inizio di ogni simulazione viene creato un file di testo in cui vengono scritte statistiche sull’oct-tree, sull’aggiornamento dei nodi, etc. Esso viene creato solo su esplicita richiesta, poiché i valori in esso memorizzati non sono necessari alla corretta interpretazione dei risultati di una simulazione. Sicuramente la sua presenza ed il suo aggiornamento peggiorano, anche se lievemente, la performance complessiva del programma.


3.5 Modifiche funzionali Sono tutti quei cambiamenti al codice atti sia a introdurre nuove funzionalità nel programma, che a togliere quelle ritenute obsolete. Le principali aggiunte sono:

• arresto forzato e riavvio di una simulazione; • supporto di file little e big endian (vedi glossario); • scelta manuale del tempo di creazione di ogni snapshot.

3.5.1 Arresto forzato di una simulazione È possibile forzare l’arresto di una simulazione in due modi:

• manualmente; • dopo un determinato tempo di CPU.

Nel primo caso è sufficiente creare un file di nome “STOP” nella directory di output del programma (specificata nel file dei parametri). Il programma, appena rileva la presenza di questo file, termina l’esecuzione scrivendo su disco, se esplicitamente richiesto, il file di restart. Questo file viene creato facendo un “dump” su disco di tutta la memoria allocata dal programma. Nel secondo caso, il programma, superato l’85 percento del tempo CPU specificato alla voce MaxAllowedCPU nel file dei parametri, inizia la procedura d’arresto della simulazione. Il 15 percento restante serve per avere la sicurezza che, se esplicitamente richiesto (si veda la voce restartFile nel file dei parametri), venga memorizzato il file di riavvio e l’eventuale snapshot finale. Se non si prevede la scrittura di un file di


riavvio la simulazione inizia la procedura d’arresto quando ha superato il 95 percento del tempo CPU.

3.5.2 Riavvio di una simulazione È possibile far ripartire una simulazione, in due modi:

• tramite il file di restart; • tramite uno snapshot.

Nel primo caso è sufficiente richiamare il programma con la seguente sintassi: inSPHector.exe -r <fileParametri> Deve essere presente nella directory di output del programma un file di restart, col nome specificato in <fileParametri>, alla voce “RestartFile”. Nel secondo caso invece si deve richiamare l’applicativo con la seguente sintassi: inSPHector.exe -r:<numeroSnapshot> <fileParametri> In questo caso deve essere presente nella directory di output del programma uno snapshot con suffisso <numeroSnapshot> e prefisso specificato in <fileParametri>, alla voce “SnapshotPrefix“.


3.5.3 Supporto di file little e big endian Poiché il programma può funzionare sia su architetture little endian che su architetture big endian, può capitare che una simulazione iniziata su un computer con una determinata architettura, debba essere ultimata su un computer con architettura diversa. Abbiamo risolto questo problema, rendendo il tutto trasparente all’utente. Il programma, infatti, riconosce automaticamente l’architettura del computer che ha generato il file delle condizioni iniziali. Se necessario i dati letti vengono convertiti nel formato corretto. È possibile inoltre specificare (nel file dei parametri alla voce “defaultSaveFormat”), se l’output del programma deve essere in formato little endian, in formato big endian o nel formato standard della macchina in cui sta girando l’applicativo. Queste funzionalità sono state attivate realizzando un modulo per l’accesso binario ai file. Per informazioni al riguardo, si rimanda all’appendice D. Prima dell’introduzione di questa possibilità, si era costretti ad usare un programma scritto in C che convertiva un file da big endian a little endian e viceversa. Inoltre, per comodità si tenevano le condizioni iniziali sia in formato big endian che in formato little endian, con un notevole spreco di risorse.


3.5.4 Scelta manuale del tempo di creazione di ogni snapshot Nel file dei parametri è necessario specificare il tempo iniziale e quello finale di una simulazione. Deve essere inoltre specificato il numero di snapshot che si desidera vengano creati durante una simulazione. Il tempo tra uno snapshot e l’altro viene calcolato, nel caso in cui la simulazione sia effettuata in coordinate fisiche, nel seguente modo:

1minmax

−−

=ratishotDesidenumeroSnap

timetimeeSnapshottempoTraDu ,

mentre se è in coordinate comoventi:

11

minmax −

=

ratishotDesidenumeroSnap

timetime

eSnapshottempoTraDu .

A volte, è necessario specificare esattamente l’istante in cui il programma deve creare uno snapshot. Per risolvere questo problema, può essere dato in input al programma un file di testo con una lista degli istanti in cui si vogliano creare gli snapshot. Il nome di questo file è specificato nel file dei parametri alla voce SnapshotTimeList.


3.6 Strutture dati e variabili principali In questa versione del programma sono state usate le “strutture” del Fortran 90. Si rimanda al glossario per maggiori dettagli. Ciò ha reso necessario una re-impostazione globale del codice il quale però è diventato molto più leggibile e manutenibile. Ogni particella è identificata da un record, di cui diamo la definizione: TYPE particle_data REAL:: Pos(0:2) ! Posizione corrente della ! particella REAL:: Vel(0:2) ! Velocità corrente della particella REAL:: Mass ! massa della particella INTEGER(INT4BYTE):: ID ! Identificatore (unico per ogni ! particella) INTEGER(INT4BYTE):: TYPE ! tipo di particella. ! 0=gas, 1=halo, 2=disk, 3=bulge, ! 4=stars REAL:: CurrentTime ! tempo corrente della particella REAL:: MaxPredTime ! massimo tempo predetto REAL:: PosPred(0:2) ! tempo predetto REAL:: VelPred(0:2) ! velocità predetta REAL:: Accel(0:2) ! accelerazione della particella REAL:: Potential ! potenziale della particella REAL:: OldAcc ! accelerazione al precedente ! timestep.

! Usato nel nuovo criterio di ! apertura.

INTEGER(INT4BYTE):: ForceFlag ! puntatore alla prossima particella ! attiva END TYPE particle_data Il programma, una volta conosciuto il numero di particelle della simulazione, alloca un array di record di tipo particle_data chiamato P, di grandezza pari al numero di particelle.


Per le particelle barioniche (le uniche che sono sottoposte anche alle leggi idrodinamiche oltre a quelle gravitazionali), è prevista un’ulteriore struttura dati. TYPE sph_particle_data REAL:: Density ! densità corrente della particella REAL:: DtDensity ! intervallo di cambimento della densità REAL:: DensityPred ! densità predetta REAL:: EgySpec ! energia interna per unità di massa REAL:: DtEgySpec ! velocità del cambiamento dell’energia ! interna REAL:: EgySpecPred ! energia interna predetta per unità di ! massa REAL:: Pressure ! pressione REAL:: Hsml ! lunghezza di smussamento REAL:: DtHsml ! velocità del cambiamento della lunghezza ! di smussamento INTEGER:: NumNgb ! numero di vicini REAL:: DivVel ! divergenza locale della velocità REAL:: CurlVel ! rotore locale della velocità END TYPE sph_particle_data

Il programma, una volta acquisito il numero di particelle barioniche della simulazione, alloca un array di record di tipo sph_particle_data chiamato SPH_P, di dimensione pari al numero di particelle barioniche. Se N è il numero di particelle barioniche della simulazione, esse devono essere presenti sia nell’array P che nell’array SPH_P. Per quanto riguarda l’algoritmo di Barnes-Hut, è stata creata una struttura di tipo nodo: TYPE Node REAL :: len ! lunghezza del lato del cubo REAL :: mass ! massa totale del nodo REAL :: s(0:2), center(0:2) ! centro di massa/

! baricentro REAL :: Q11, Q22, Q33, Q12, Q13, Q23, P! tensori del quadrupolo REAL :: vs(0:2) ! velocità del centro di massa


REAL :: tilu ! tempo di ultima mod. del nodo REAL :: oc ! variabile per il nuovo

! criterio di apertura REAL :: hmax ! massima lunghezza di

! smussamento INTEGER(int4byte) :: suns(0:7) ! discendenti INTEGER(int4byte) :: sibling ! contiene il prossimo nodo, se ! questo nodo non deve essere

! ulteriormente aperto INTEGER(int4byte) :: partind ! indice della prima ! particella nel nodo INTEGER(int4byte) :: count ! numero totale di particelle ! di un nodo END TYPE Node

L’oct-tree viene costruito a partire da un array di nome nodes con elementi di tipo Node

Servono degli array di tipo “integer”: next: serve a trovare tutte le particelle appartenenti ad un nodo. nextnode: serve a trovare il prossimo nodo dell’albero. father: contiene l’indice nell’array del genitore di quel nodo. Per quanto riguarda l’integrazione temporale un array di interi viene interpretato come un albero binario di ricerca. Per maggiori dettagli si rimanda al paragrafo 3.11.2.


3.7 Implementazione dell’algoritmo di Barnes-Hut

3.7.1 Introduzione L’algoritmo di Barnes-Hut è disponibile in rete in molte versioni. Una versione scritta in C prevede l’allocazione dinamica della memoria e dei nodi che compongono l’oct-tree. Purtroppo questa versione presenta due svantaggi:

• l’allocazione dei nodi richiede tempo macchina. Poiché durante una simulazione, l’albero viene ricostruito più volte, l’overhead portato dalla allocazione dinamica dei nodi può non essere trascurabile;

• ogni nodo può avere come figli sia un altro nodo che una particella. Ciò è possibile in C attraverso l’uso di un typecasting esplicito, rendendo però il codice meno leggibile.

Nel nostro codice abbiamo usato un array per contenere i nodi dell’albero, eliminando gli inconvenienti sopraccitati. Non potendo sapere a priori il numero di nodi che servono per costruire l’albero, il dimensionamento di questo array è problematico. V.Springel, N.Yoshida e S.D.M.White (2000), compiendo numerose prove, mostrano che, se N è il numero di particelle in una particolare simulazione, generalmente bastano N*0.7 nodi. Abbiamo utilizzato questa stima nel nostro codice. Quando si supera questo valore, viene effettuata una re-allocazione dinamica dell’array contenenti i nodi. Per maggiori dettagli si veda il paragrafo 3.7.4.


3.7.2 Implementazione La nostra implementazione dell’albero di Barnes-Hut prevede l’uso di array. Per capire se un figlio di un nodo è ancora un nodo o una particella, abbiamo usato un piccolo accorgimento. Sia numPart il numero di particelle del sistema. Poiché l’array delle particelle deve avere lunghezza numPart e ogni particella ha un indice che varia da 1 a numPart, abbiamo assegnato indice numPart al primo elemento dell’array dei nodi. In pseudo-codice abbiamo: MaxNodes=0.7*NumPart alloca array_delle_particelle(1..NumPart) alloca array_dei_nodi (NumPart+1:NumPart+MaxNodes) Per indicare l’assenza di figli di un nodo abbiamo utilizzato il valore -1. L’algoritmo di Barnes-Hut prevede che ogni nodo o foglia dell’albero debba conoscere il genitore. Per fare questo abbiamo viene allocato un array di dimensione NumPart+MaxNodes, che memorizza questa informazione. La nostra implementazione dell’algoritmo di Barnes-Hut, prevede la costruzione di un albero per ogni tipo di particella della simulazione. Questo velocizza notevolmente il codice, poiché ad un certo passo temporale solo una piccola percentuale di particelle si muove. Inoltre tipi di particelle diverse si possono muovere con tempi diversi. Supponiamo ad esempio di far partire una simulazione con N particelle barioniche, e 10N particelle di materia oscura. Supponiamo, che ad un certo passo temporale, solamente N/10 particelle barioniche debbano essere mosse.


Figura 6: esempio di multi-albero e relativa rappresentazione con array

In questo caso conviene ricostruire solamente un albero contenente come foglie N particelle, piuttosto che uno contenente come foglie 10N+N particelle. Inoltre molte volte succede che due particelle di tipo diverso hanno posizione spaziale quasi uguale In tal caso un algoritmo che prevede un solo albero, necessita di una profondità notevole, con ripercussioni sul tempo impiegato nella visita.


3.7.3 Aggiornamenti dinamici dell’albero Se il numero di particelle che devono essere fatte avanzare ad un certo timestep è piccolo, la predizione della posizione di tutte le particelle e la ricostruzione dell’albero possono portare ad un significativo overhead. D’altra parte il modo in cui le particelle sono raggruppate in una gerarchia evolve in maniera relativamente lenta nel tempo. È perciò sufficiente ricostruire tale raggruppamento solo ogni tanto, a patto che si possano ottenere accurate informazioni sul nodo (centro di massa, momenti multipolo, etc.) nell’istante attuale. Nel codice si ottiene questo risultato predicendo “al volo” le proprietà dell’albero, invece di predire la posizione di tutte le particelle ad ogni timestep. Ciascun nodo contiene una velocità del centro di massa, in aggiunta alla sua posizione all’istante in cui è stato costruito. La nuova posizione dei nodi può essere predetta mentre l’albero viene esplorato e solo le posizioni dei nodi che sono visitati viene predetta. Nel nostro algoritmo trascuriamo una possibile variazione nel tempo del momento di quadrupolo dei nodi. Abbiamo comunque introdotto un meccanismo addizionale che reagisce a veloci variazioni dei nodi dell’albero. Nel momento in cui il centro di massa di un nodo si sposta più di una certa frazione del lato del nodo stesso, il nodo viene completamente aggiornato, ossia il centro di massa, la velocità del centro di massa ed il momento quadrupolo vengono ricalcolati. Verrà tarata la lunghezza del lato dei nodi dell’albero se qualcuna delle sue particelle avrà perso il suo nodo d’appartenenza originario. Tipicamente si ricostruisce l’albero quando sono stati effettuati un totale di circa N1.0 calcoli di forze dalla sua ultima piena ricostruzione.


Nel file dei parametri è possibile scegliere questo valore alla voce TreeUpdatingLimit. Se si sceglie come valore 0, l’albero viene sempre ricostruito. Con questo criterio la costruzione dell’albero è solo una frazione non significativa del tempo di calcolo totale, mentre non è stata riscontrata una significativa perdita di accuratezza delle forze.

3.7.3.1 Complessità Poiché ad un certo timestep solo un piccolo numero M di particelle viene fatto avanzare, il tempo necessario per l’aggiornamento dell’albero è

)lg( 8 NMO .

3.7.4 Riallocazione dinamica dei nodi Come già accennato, non possiamo conoscere a priori quanta memoria dobbiamo allocare per l’array che rappresenta i nodi dell’oct-tree. Tuttavia generalmente un numero di nodi pari al settanta percento del numero di particelle, basta per una simulazione. Può succedere, però, che l’array sia troppo piccolo. In tal caso, il programma crea un altro array supplementare che diventerà il nuovo array dei nodi, il dieci percento più grande di quello che contiene l’oct-tree. Successivamente copia tutti i valori dal vecchio array a quello nuovo e infine disalloca quello vecchio.


3.8 Nuovo criterio per l’apertura dell’oct-tree L’accuratezza delle forze risultanti percorrendo l’albero, dipende sensibilmente dal criterio utilizzato per decidere se l’approssimazione multipolo per un certo nodo è accettabile o se il nodo deve essere “aperto” per un ulteriore raffinamento. Come già detto, ogni nodo rappresenta un regione cubica nello spazio di simulazione in cui vi sono una o più particelle. Se l’approssimazione è accettabile, allora si dice che le particelle del nodo sono “ben separate” rispetto alla particella su cui si vuole valutare la forza. L’apertura standard, ovvero quella proposta da Barnes-Hut, si basa su un semplice criterio geometrico. L’idea è quella di valutare la proporzione tra il lato l di un nodo e la distanza r tra il centro di massa delle particelle incluse nel nodo e la particella su cui si vuole valutare la forza. Se:

Θ<rl

(3.2)

dove Θ è una costante fissata a priori e minore di uno, le particelle sono “ben separate”.


l = lato del cubo

c.m.

r =distanza tra la particella ed ilcentro di massa del cubo

particella su cuicalcolare la forza

Figura 7: approssimazione di Burnes-Hut

L’apertura standard, ovvero quella proposta da Barnes-Hut, stima l’errore ad ogni iterazione nodo-particella, confrontando il termine quadrupolo,

4

2

rml

, con il termine monopolo, 2r

m.

Come riportato da Salmon & Warren (1994), nel caso peggiore il comportamento di questo criterio di apertura è poco affidabile. Sebbene accada di raro nelle simulazioni astrofisiche, il criterio di Barnes-Hut può occasionalmente portare ad errori di approssimazione molto grandi nel calcolo delle forze. Un altro problema si presenta quando lo si utilizza con un alto redshift. In questo caso il campo di densità è prossimo all’omogeneità e l’accelerazione di ogni particella è piccola. Quindi la forza risultante è calcolata tramite un processo di cancellazione tra forze parziali relativamente grandi. Anche un piccolo errore relativo, in questo caso, è sufficiente per provocare un grande errore nel calcolo della forza risultante.


Per una distribuzione omogenea delle particelle, il criterio di apertura di Barnes-Hut richiede un piccolo valore dell’angolo di apertura. In una simulazione cosmologica, il valore delle forze tra una particella data e l’albero, per una certa dimensione dell’angolo di apertura può variare di parecchi ordini di grandezza. Come suggerito da Salmon & Warren (1994), è conveniente usare un criterio di apertura che limiti l’errore assoluto in ogni iterazione cella-particella. Il nostro approccio prevede che si conosca la dimensione della forza assoluta prima di iniziare il cammino dell’albero. Nel nostro codice usiamo l’accelerazione del precedente timestep come comodo valore approssimato. Imponiamo che la stima dell’errore per una accettabile approssimazione multipolo sia una qualche piccola frazione, α , di questa forza totale. Secondo V.Springel, N.Yoshida e S.D.M.White (2000) ne segue che un nodo è aperto se:

64 raMl oldα>

Studi accurati hanno mostrato che questo criterio è più efficiente del criterio di Barnes-Hut poiché, a parità di costo computazionale, produce un calcolo delle forze più accurato. Inoltre lo stesso valore α produce una confrontabile accuratezza nel calcolo delle forze indipendentemente dalla distribuzione delle particelle.


3.9 Supporto delle coordinate comoventi Poiché l’universo si espande, la sua dimensione tipica, a, aumenta nel tempo. Anche la distanza r tra due galassie lontane aumenta nel tempo in modo proporzionale ad a: ar ∝ . Per studiare l’evoluzione dinamica delle strutture cosmiche, è utile mettersi in un sistema di riferimento che “rimuova” la pura espansione geometrica dell’universo. Come esempio, pensiamo ad un palloncino di gomma su cui sono disegnati meridiani e paralleli. Gonfiando il palloncino, anche il sistema di meridiani e paralleli si espande. Un sistema di riferimento basato su meridiani e paralleli viene chiamato comovente. La distanza comovente x tra due galassie lontane è legata a quella fisica (o propria) r dalla semplice relazione xar ⋅= . In un sistema di riferimento comovente il moto relativo delle galassie è determinato unicamente dalla mutua attrazione gravitazionale. Tuttavia, è nel sistema di riferimento proprio o fisico che l’energia di un sistema isolato (come ad esempio l’universo) si conserva. Di conseguenza, le equazioni del moto nel sistema comovente andranno modificate per tener conto dell’espansione dell’universo. L’equazione che descrive l’espansione dell’universo è quella di Friedmann per il fattore di scala:

( ) ( )3

18 110

22

−− +−Ω=

=

aGdtda

ta iρπ (3.3)

in cui t è il tempo proprio, a(t) è una lunghezza da determinare, detta fattore di scala o parametro di espansione, G è la costante di gravitazione, 0ρ è la

densità al tempo attuale. Ω è la percentuale di materia ed energia che esiste


nell’universo in rapporto alla massa critica 2)8

3( H

Gc πρ = , dove H è la

costante di Hubble (circa 565 ± (km/s)/Mpc). Per quanto riguarda le equazioni di Eulero, se agisce solo la gravità, posto x=r/a (che definisce la coordinata comovente x), si ricava:

i

i

jiij

ijjii m

aGm

aa F

x

xv

av 3

332 −

≠

− =−== ∑ (3.4)

Per ulteriori dettagli si rimanda a [Dal Borgo, 2000].

3.10 Supporto delle condizioni periodiche al contorno Compilando il codice con l’opzione di precompilazione PERIODIC, vengono supportate le condizioni periodiche al contorno. In questo caso abbiamo un calo delle prestazioni del programma che viene appesantito da numerosi calcoli. Se il dominio di una simulazione è un cubo di lato L, le condizioni periodiche al contorno implicano che:

• se una particella esce da un lato del cubo della simulazione, essa deve rientrare dal lato opposto;

• una particella, oltre all’accelerazione dovuta a tutte le altre particelle della simulazione, subisce anche l’influenza di tutti i corpi presenti nel resto dell’universo.

Per implementare questa caratteristica, nel codice abbiamo fatto uso dei metodi di Ewald. Per una loro dettagliata descrizione si rimanda a [Dal Borgo, 2000].


3.11 Integrazione temporale Come integratore temporale, abbiamo utilizzato una variante del metodo Leap-Frog con passo esplicito. Iniziamo con il descrivere l’integratore temporale per una singola particella.

3.11.1 Integratore temporale per una singola particella La posizione di una particella a metà di un timestep t∆ è predetta con:

2~ )()()

21

( tnnn ∆+=

+vrr (3.5)

e l’accelerazione in questa posizione diventa:

)21(~

)21

(n

+Φ−∇=

+

nr

a (3.6)

La particella avanza in accordo alla relazione:

tnnn ∆+=

++ )21

()()1( avv , (3.7)

[ ] tnnnn ∆++= ++ )1()()()1(

21

vvrr . (3.8)

In questo schema l’intervallo temporale può variare da passo a passo.


Nel codice è possibile scegliere di calcolare t∆ in due modi diversi, a seconda che il valore del campo TimestepType nel file di configurazione sia zero o uno. Se TimestepType è uguale a zero abbiamo:

aεα '2 tolt =∆ (3.9)

mentre nel secondo:

a

''tolt

α=∆ (3.10)

dove: 'tolα e ''

tolα sono parametri adimensionali di accuratezza e vengono letti dal file di configurazione rispettivamente alla voce DeltatTollerance1 e DeltatTollerance2; ε è la lunghezza di smussamento gravitazionale e a è l’accelerazione della particella presa in esame. Secondo V.Springel, N.Yoshida e S.D.M.White (2000) la (3.10) dà risultati migliori di molte altre scelte, compresa la (3.9), rimasta nel codice perché essa costituisce l’approccio classico. Per quanto riguarda le particelle SPH, poiché esse sono soggette anche a forze di tipo idrodinamico, si devono leggermente modificare le equazioni per il calcolo dei passi temporali. Sia i una particella:

( ) )6.01)(,max( visciiii

icouri ch

ht

αα

++⋅∇=∆

vv (3.11)

dove viscα è un parametro di viscosità artificiale, courα è un parametro di

accuratezza, chiamato fattore di Courant.


3.11.2 Integratore temporale per sistemi N-Body Nel codice abbiamo implementato un integratore a timestep completamente flessibili, simile a quello impiegato da Groom (1997) e Hiotelis & Voglis(1991). Con questo schema lo stato dinamico di una particella può essere predetta al tempo [ ]ii ttt ∆±∈ 5.0 con accuratezza del primo ordine. Ne deriva, quindi, che la particella da muovere ad un certo istante t è quella che ha il minimo tempo predetto )5.0min( iip tt ∆+≡τ .

Il tempo pτ diventa il nuovo tempo corrente della simulazione.

Per far avanzare la corrispondente particella per prima cosa prediciamo le posizioni di tutte le particelle al tempo pτ in accordo con:

)(~

ipiii t−+= τvrr (3.12)

Basandoci su queste posizioni, l’accelerazione di una particella di indice k nel mezzo di questo timestep è calcolata con:

kri

n

k ~)

21

()~(ra Φ−∇=

+ (3.13)

La posizione e la velocità della particella di indice k sono rispettivamente:

[ ] )()1()()()1(kp

nk

nk

nk

nk t−++= ++ τvvrr (3.14)

dove

)(2)

21

()()1(kp

n

kn

kn

k t−+=++ τavv (3.15)


Il suo tempo corrente diventa:

(2+= knuovok tt )kp t−τ (3.16)

Naturalmente è poco pratico avanzare una singola particella ad ogni nuovo istante, perché il calcolo e gli aggiornamenti dinamici dell’albero introducono un overhead. Per questa ragione avanziamo più particelle contemporaneamente. Le particelle vengono ordinate rispetto ai loro istanti di predizione:

iipred

i ttt ∆+=21

. (3.17)

La simulazione avanza sempre la particella con il più piccolo predit e tutte

quelle che soddisfano:

iip tt ∆+≥41

τ . (3.18)

Questa condizione seleziona un gruppo di particelle al limite inferiore della linea temporale e a tutte le particelle del gruppo è garantito l’avanzamento di al più metà del massimo timestep permesso loro. L’esperienza dimostra che il gruppo che è avanzato ad un determinato passo, spesso contiene solo una piccola frazione del numero totale di particelle. E necessario eliminare ogni overhead dovuto al calcolo temporale. Abbiamo sicuramente bisogno della particella con il minimo istante predetto e, inoltre, di tutte le particelle che la seguono nella linea temporale. Questo problema è stato risolto, mantenendo i tempi massimi predetti in un albero rosso-nero. Ricordiamo che un albero rosso-nero è un albero binario di ricerca, approssimativamente bilanciato. Per maggiori dettagli si rimanda all’appendice C. In questo modo, la ricerca della particella con il minimo istante predetto costa )(log NO . Inoltre le particelle che sono state fatte avanzare possono essere tolte e reinserite nell’albero sempre con un costo )(log NO .


3.12 Ricerca dei vicini

Nell’algoritmo SPH un’operazione base consiste nel cercare i “vicini” di ciascuna particella, entro un certo raggio ih .

Analogamente alla gravità l’algoritmo più semplice, ovvero quello che controlla la distanza di tutte le coppie di particelle, ha ordine )( 2NO . Ciò rallenta notevolmente l’elaborazione nel caso in cui vi siano molte particelle. Se la distribuzione delle particelle è omogenea, l’algoritmo più efficiente è quello che si basa su di una griglia di ricerca, in cui ogni cella ha dimensione pari al raggio di ricerca. Tramite delle liste si potrebbe, poi, arrivare solo alle particelle che giacciono in una specifica cella della griglia e con una semplice ricerca ad intervalli trovare le particelle che stanno entro il raggio r. Purtroppo però, per distribuzioni altamente non omogenee e variando il raggio ih , il precedente approccio degrada velocemente poiché il reticolo

scelto per la griglia non ha la migliore dimensione per tutte le particelle. Un approccio alternativo e flessibile prevede un albero di ricerca geometrico. Fortunatamente l’oct-tree si presta molto bene a questo scopo. La ricerca dei vicini è effettuata visitando l’albero. Una cella è aperta se ha una sovrapposizione spaziale con un intervallo di ricerca rettangolare. Se si arriva ad una cella con una sola particella questa è aggiunta alla lista di iterazione solo se sta dentro al raggio di ricerca. Nell’algoritmo appena descritto, l’albero necessita di una esplorazione che raggiunga le foglie contenenti le particelle che andranno mosse. Abbiamo


bisogno, durante la costruzione dell’albero, di creare una lista che ci permetta un reperimento di tutte le particelle appartenenti ad una data cella. Nel nostro programma, seguendo il lavoro di V.Springel, N.Yoshida e S.D.M. White (2000), abbiamo usato uno stratagemma aggiuntivo. Un cammino sull’albero viene terminato se la cella giace completamente all’interno dell’intervallo di ricerca. In tal caso tutte le particelle appartenenti a quella cella possono essere aggiunte alla lista di quelle che debbono essere mosse. Con una leggera modifica dell’algoritmo di attraversamento dell’albero possiamo, inoltre, individuare tutte le particelle con ),max( jiij hhr < .

A tale scopo memorizziamo in ciascun nodo dell’albero la massima lunghezza di smussamento nodoh che si presenta tra queste particelle. La verifica per la sovrapposizione è quindi semplicemente eseguita tra un cubo di dimensione ),max( nodoi hh centrato nella particella i ed il nodo

stesso, dove nodoh è la massima lunghezza di smussamento tra tutte le

particelle del nodo. Rimane quindi da rendere il numero di vicini attorno ad una data particella SPH approssimativamente (o esattamente) costante. Ciò si risolve predicendo un valore ih

~ per la lunghezza di smussamento.

Il suo calcolo viene basato sul precedente timestep; viene utilizzata la lunghezza di smussamento del precedente timestep, il numero di vicini iN

in quel timestep e la divergenza della velocità locale:

thNN

hh ii

svecchioii ∆+

+= &3

1

)( 121~

.


dove iii hh )*(31

v∇=& e sN è il numero di vicini desiderato.

Il termine nelle parentesi quadre è volto a riportare il numero di vicini al valore originario desiderato nel caso in cui iN se ne allontani.

Nel caso in cui, nonostante ciò, il numero di vicini cada al di fuori di un predefinito intervallo di tolleranza, iterativamente ih viene modificato al

fine di riportare il numero di vicini all’intervallo desiderato. Il codice permette, tramite un’apposita opzione nel file di configurazione alla voce MinSPHtoGravSoftRatio, di impostare una minima lunghezza di smussamento per le particelle barioniche. Vengono ammessi un numero di vicini “smussati” maggiore di sN se si

raggiunge questo valore minimo. Lo schema seguito dal programma è il seguente. Si procede, come primo passo, alla ricerca dei vicini entro un raggio

ihR 2.1= .

Di questi vicini selezioniamo un insieme con cardinalità sN .

Se ci sono meno di sN particelle nell’intervallo di ricerca, o se la distanza

della particella più vicina è più grande di R, la ricerca viene ripetuta per un più ampio intervallo di ricerca. Al primo timestep non si conosce il precedente ih , quindi noi seguiamo l’albero dei vicini all’indietro, ovvero dalla particella che consideriamo, fino a quando riusciamo ad ottenere una prima ragionevole stima della densità locale delle particelle (basato sul numero N di particelle in un nodo di volume 3l ). Questo fornisce una stima iniziale per ih .


3.13 Moduli utilizzati Come di può vedere dalla figura 8, il programma è costituito da cinque moduli principali:

• I/O su file; • Variabili, costanti e funzioni globali; • SPH: si occupa della fluido-dinamica; • Algoritmo di Barnes-Hut: gestisce l’albero di Barnes-Hut, calcola la

forza ed il potenziale di ogni particella; • Integratore temporale: è la parte del programma che si incarica di

far avanzare le particelle. Begrun_m si occupa dell’inizializzazione della simulazione, mentre run_m la esegue. In main c’è l’entry point del programma. Il programma inizia dal file main.f90. Viene analizzata la riga di comando e viene visualizzato a video ogni errore commesso dall’utente. Se l’utente desidera avviare (o riavviare) una simulazione, viene dapprima chiamato il metodo begrun, che si incarica di preparare l’ambiente e poi il metodo run, che avvia effettivamente la simulazione. I moduli inclusi sono: begrun_m, run_m.


Figura 8: schema dei moduli e delle relative inclusioni


Di seguito elenchiamo, in ordine alfabetico, i moduli del programma e una loro sintetica descrizione. Inoltre verranno elencati i corrispettivi moduli inclusi, tutte le routine, le variabili, le costanti e le interfacce pubbliche. Modulo: begrun_m.f90 Descrizione: In questo modulo sono presenti funzioni utili all’avvio

corretto di una simulazione, come la lettura del file di configurazione, l’apertura dei file di output, etc.

Moduli utilizzati: IO_readic_m, IO_restart_m, TI_timestep_m, TI_accel_m Subroutine: begrun Modulo: binary_m.f90 Descrizione: in questo modulo sono presenti funzioni e subroutine per

l’accesso binario e sequenziale ai file, più altre funzioni d’utilità.

Moduli utilizzati: nessuno (il modulo deve essere utilizzabile anche in altri contesti) Interfacce: binWrite, binRead Funzioni: machineArch, typeOfNotFormattedFile, fileLength Subroutine: binOpen, binClose, binFlush, binRewind, binConversion Modulo: BH_force_m.f90 Descrizione: modulo per il calcolo delle forze gravitazionali. Moduli utilizzati: BH_tree_m, TI_predict_m Subroutine: compute_gravity_forces


Modulo: BH_periodic_m.f90 Descrizione: modulo per il calcolo delle condizioni periodiche al contorno Moduli utilizzati: GLOBAL_m, erf_m, binary_m Funzioni: ewald_pot_corr Subroutine: ewald_init, ewald_force, ewald_corr Modulo: BH_potential_m.f90 Descrizione: Modulo per il calcolo del potenziale di ogni particella Moduli utilizzati: BH_tree_m, TI_predict_m Subroutine: compute_potential Modulo: BH_tree_m.f90 Descrizione: in questo modulo sono presenti le funzioni per la gestione

dell’oct-tree e del relativo profiling. Moduli utilizzati: BH_treevar_m, BH_periodic_m Funzioni: BH_tree_build, BH_get_cost_part_part_inter,

BH_get_cost_part_node_inter Subroutine: BH_allocate, BH_reset_cost, BH_ngb_update_nodes,

BH_update_node Modulo: BH_treevar_m.f90 Descrizione: in questo modulo sono presenti tutte le variabili necessarie

alla creazione ed alla gestione dell’oct-tree. Poiché ci sono solo variabili, esse devono essere tutte pubbliche, altrimenti la loro presenza non sarebbe necessaria.

Moduli utilizzati: GLOBAL_m Nuovi tipi: Node Costanti: KERN_LEN=10000


Variabili: nodes, nodes_new, next, nextnode, nextnode_new, father, father_new, Part, SphPart , S_of_a_inverse, BoxHalf, Box, trees, numnodestree, numnodestotal , MaxNodes, treecost, treecost_quadru, knlrad, knlforce, knlpot, knlW2, knlW3, knlW4, last

Modulo: erf_m.f90 Descrizione: in questo modulo sono presente funzioni che restituiscono

l’integrale della funzione errore ed il suo complementare. Funzioni: erf, erfc Modulo: GLOBAL_const.f90 Descrizione: in questo modulo sono presenti tutte le costanti necessarie

alla simulazione. Costanti: DOUBLE=8,INT4BYTE=4, MAX_REAL_NUMBER=1e37

MIN_REAL_NUMBER = 1e-37, THIRD = 1.0 /3.0 PI = 3.14159265358979323846, PI_INV = 1.0/PI, SECFACTOR = 1.2_DOUBLE, LN2 = 0.69314718 GRAVITY = 6.672e-8, SOLAR_MASS = 1.989e33, SOLAR_LUM = 3.826e33, RAD_CONST = 7.565e-15, AVOGADRO = 6.0222e23, BOLTZMANN = 1.3806e-16, GAS_CONST = 8.31425e7, C = 2.9979e10, PLANCK = 6.6262e-27, CM_PER_MPC = 3.085678e24, PROTONMASS = 1.6726e-24, HUBBLE = 3.2407789e-18, SEC_PER_MEGAYEAR = 3.155e13, SEC_PER_YEAR = 3.155e7, GAMMA = 5.0_DOUBLE/3.0_DOUBLE, GAMMA_MINUS1 = 2.0_DOUBLE/3.0_DOUBLE, KERNEL_TABLE = 1000, MAX_NGB = 20000, MAXLEN_SNAPSHOTTIMELIST = 350, DEFAULTTREEALLOCFACTOR =.7 FdInfo=100, FdEnergy=102, FdTimings=104, FdCPU=106


Modulo: GLOBAL_m.f90 Descrizione: in questo modulo sono presenti subroutine generali necessarie

al programma. Moduli utilizzati: GLOBAL_set_m, various_m Subroutine: allocate_memory, free_memory, find_next_outputtime,

compute_global_quant_of_system, printPrecompiledOptions Modulo: GLOBAL_set_m.f90 Descrizione: Funzioni per l’inizializzazione del kernel, delle lunghezze di

smussamento, etc. Moduli utilizzati: GLOBAL_vars_m Subroutine: set_units, set_sph_kernel, set_softenings Modulo: GLOBAL_vars_m.f90 Descrizione: in questo modulo sono presenti tutte le variabili necessarie

alla simulazione. Poiché non ci sono né subroutine né funzioni, le variabili sono tutte pubbliche, altrimenti la loro presenza non sarebbe necessaria.

Moduli utilizzati: GLOBAL_const_m Nuovi tipi: global_data_all_processes, particle_data, sph_particle_data, timetree_data, state_of_system, io_header Variabili: NumForceUpdate, NumSphUpdate, IndFirstUpdate,

Num_nodeupdates, Num_nodeupdate_particles, RestartFlag, RestartSnapshotFlag,EnergyRestart, NumPart, N_gas, ParameterFile, Kernel, KernelDer, KernelRad, Simulation, P, SPH_P, PTimeTree, SysState, SysStateAtStart, SysStateAtEnd, CPUThisRun, header, snapversion, machineArchitecture


Modulo: io_m.f90 Descrizione: modulo che si incarica della lettura e della scrittura delle

statistiche, del file dei tempi degli snapshot, delle intestazioni dei file, etc.

Moduli utilizzati: GLOBAL_m, binary_m Subroutine: open_outputfiles, close_outputfiles, read_outputlist,

infoInfoFile, infoCpuFile, infoTimingsFile, infoEnergyFile, writeSimulationVar, readSimulationVar, readParticle_data, writeParticle_data, readSphParticleData, writeSphParticleData, readHeader, writeHeader

Modulo: IO_readic_m.f90 Descrizione: Contiene delle routine per la lettura dei file delle condizioni

iniziali. Moduli utilizzati: IO_m Subroutine: read_ic_type_1, read_ic_type_3, read_ic_type_4,

Read_ic_type_5 Modulo: IO_restart_m.f90 Descrizione: Contiene le routine per il restart di una simulazione. Moduli utilizzati: IO_m Subroutine: fullRestart, snapshotRestart, saveSnapshot


Modulo: run_m.f90 Descrizione: Contiene la funzione run, incaricata di eseguire l’intera

simulazione Moduli utilizzati: BH_potential_m, IO_restart_m, TI_timestep_m,

TI_accel_m Subroutine: run Modulo: SPH_cooling_m.f90 Descrizione: in questo modulo sono presenti delle funzioni per il calcolo

del cooling. Moduli utilizzati: GLOBAL_m Funzioni: DoCooling, AbundanceRatios Subroutine: InitCool, IonizeParams Modulo: SPH_density_m.f90 Descrizione: in questo modulo sono presenti delle funzioni per il calcolo

della densità di una particella barionica. Moduli utilizzati: SPH_neighbour_m Subroutine: density Modulo: SPH_hydra_m.f90 Descrizione: modulo per il calcolo delle forze idrodinamiche Moduli utilizzati: SPH_neighbour_m, SPH_cooling_m (incluso se si

compila con la direttiva COOLING) Subroutine: SPH_compute_forces


Modulo: SPH_neighbour_m.f90 Descrizione: modulo per il calcolo dei vicini delle particelle barioniche Moduli utilizzati: BH_tree_m Funzioni: ngb_treefind, ngb_treefind_variable, ngb_treefind_pairs Subroutine: ngb_treeallocate, ngb_setup_smoothinglengths Modulo: TI_accel_m.f90 Descrizione: in questo modulo è presente una funzione che calcola

l’accelerazione per tutte le particelle attive. Moduli utilizzati: SPH_density_m, SPH_hydra_m, BH_force_m Subroutine: compute_accelerations Modulo: TI_predict_m.f90 Descrizione: modulo per l’avanzamento delle particelle Moduli utilizzati: GLOBAL_m Subroutine: predict, predict_sph_particles, advance,

predict_collisionless_only Modulo: TI_timeline_m.f90 Descrizione: modulo per il calcolo delle particelle da spostare ad un certo

timestep Moduli utilizzati: GLOBAL_m Variabili: root Subroutine: find_next_time, construct_timetree, timetree_insert_node,

timetree_delete_node


Modulo: TI_timestep_m.f90 Descrizione: modulo che si incarica di determinare i timestep. Moduli utilizzati: TI_timeline_m Subroutine: find_timesteps Modulo: various_m.f90 Descrizione: modulo ad utilizzo generale che contiene routine di utilità Interfacce: swap Funzioni: parse, uppercase, lowercase, second, timediff, writeString Subroutine: endrun


3.15 Valutazione dell’efficienza

In questo paragrafo presentiamo i risultati di due simulazioni eseguite sia con il codice originale [D. Dal Borgo, 2000] e con quello reimpostato, partendo da identiche condizioni iniziali. Non abbiamo potuto confrontare i codici usando un sistema di riferimento comovente, visto che tale possibilità non è prevista nella versione del programma da cui siamo partiti. Le nostre simulazioni evolvono sistemi chiusi, la cui energia totale rimane costante nel tempo. Gli algoritmi da noi utilizzati, però, compiono delle approssimazioni numeriche; un programma è tanto più preciso quanto meglio conserva l’energia, ovvero quanto più piccolo è l’errore relativo ε :

ε=∆

=−

ii

i

EE

EEE

dove iE è l’energia iniziale del sistema ed E è quella corrente.

Sia origε l’errore ottenuto con il codice originale ed reimpε l’errore ottenuto

con il codice reimpostato. Sia origT il tempo di CPU speso dal programma originale e reimpT il tempo di

CPU speso dal programma reimpostato. Siano finaleorig,ε e finalereimp,ε gli errori al termine di una simulazione, ottenuti

rispettivamente con il codice originale ed il codice reimpostato. Analogamente siano finaleorigT , e finalereimpT , i tempi di CPU spesi

rispettivamente dal codice originale e dal codice reimpostato.


Figura 9: Simulazione 1: tempi di calcolo.

finaleorigfinalereimp ,, 9.0 εε = (ad ogni passo origreimp εε < ).

La macchina su cui abbiamo effettuato i test è basata su un processore AMD Athlon 650MHz, con 128MByte di RAM e sistema operativo Windows 2000. Il compilatore utilizzato è il DIGITAL Visual Fortran versione 5.0d. La prima simulazione, 1S , utilizza quattromila particelle delle quali duemila sono di gas e duemila di materia oscura. Il numero dei vicini nell’algoritmo SPH è impostato a 40.


Per studiare il comportamento del codice reimpostato all’aumentare del numero di particelle e dei vicini, abbiamo eseguito una seconda simulazione, 2S . Essa utilizza seimila particelle, delle quali tremila sono di gas e tremila di materia oscura. Il numero dei vicini nell’algoritmo SPH è impostato a 80. Le condizioni iniziali prevedono che le particelle siano distribuite uniformemente all’interno di una sfera. Fissati i dati iniziali, abbiamo suddiviso l’intervallo temporale di ogni simulazione in 50=P passi. Sia per 1S che per 2S abbiamo fatto evolvere il sistema con il programma originale eseguendo un “run” R. Con quello reimpostato abbiamo effettuato due run. Nel primo, εR , abbiamo

deteminato i parametri del codice in modo che risulti finaleorigfinalereimp ,, εε ≅ .

Nel secondo run, TR , abbiamo determinato i parametri in modo che

risultasse finaleorigfinalereimp TT ,, ≅ .

La figura 9 si riferisce alla prima simulazione 1S . Essa confronta i risultati di R ed εR .

Esaminando la figura si vede che per effettuare i primi 35 passi della simulazione il tempo di CPU è pressoché lo stesso con i due programmi. Dopo il trentacinquesimo passo, reimporig TT >> ; origT è fino a 8 volte più

grande di reimpT . Il forte aumento di origT per 35>P corrisponde a forti

accelerazioni nel sistema simulato. In presenza di queste accelerazioni, i passi temporali diventano molto piccoli e conseguentemente vi è un aumento del tempo di calcolo necessario per eseguire ogni passo. La figura 10, sempre relativa a 1S , confronta i risultati di R ed TR .

In questo caso abbiamo finaleorigfinalereimp TT ,, 8.0= . Abbiamo riscontrato

origreimp TT < ad ogni passo.


Figura 10: Simulazione 1: conservazione dell’energia.

finaleorigfinalereimp TT ,, 8.0= (ad ogni passo origreimp TT < ). Dalla figura 10 si vede che il codice reimpostato è fino a 7 volte più preciso del codice da cui siamo partiti. Il forte aumento di origT per 35>P corrisponde, come già osservato, a forti

accelerazioni nel sistema simulato. In presenza di queste accelerazioni, il nuovo criterio di apertura dell’oct-tree nel codice reimpostato, permette di ottenere una precisione nettamente migliore rispetto all’apertura standard di Barnes-Hut, implementata nel codice originale. Quindi, dalle figure 9 e 10 si evince che in queste simulazioni il codice reimpostato gira fino a 8 volte più rapidamente se si impone


finaleorigfinalereimp ,, εε ≅ , ed è fino a 7 volte più accurato se si impone

finaleorigfinalereimp TT ,, ≅ .

Figura 11: Simulazione 2: tempi di calcolo.

finaleorigfinalereimp ,, 85.0 εε = (ad ogni passo origreimp εε < ).

Per verificare che il codice reimpostato sia più efficiente anche quando aumenta il numero di particelle ed il numero di “vicini”, nella seconda simulazione, 2S , abbiamo aumentato i due valori. La figura 11 si riferisce a 2S e confronta R con εR . In questa simulazione

reimpT è fino a 7 volte inferiore a origT .


Figura 12: Simulazione 2: conservazione dell’energia

finaleorigfinalereimp TT ,, 9.0= (ad ogni passo origreimp TT < ). La figura 12 mostra che il codice in 2S , reimpε in TR è fino a 6 volte

inferiore a origε .

Si noti che, ad ogni passo, anche in questo caso, in εR abbiamo

origreimp TT < e in TR , origreimp εε < .

Questi risultati confermano che il codice reimpostato è più efficiente dell’originale anche quando si aumentano il numero delle particelle e dei vicini.

CAPITOLO 4. UN’UTILITY: STOB 77

4. L’utility StoB

4.1 Introduzione L’implementazione di inSPHector è stata lunga e faticosa (ci si è trovati davanti a numerosi problemi, tra cui alcuni legati al fatto che il codice da cui siamo partiti non era omogeneo, scritto da molte persone e in tempi diversi). Ecco che a volte modifiche anche minimali, portavano ad errori nei risultati anche molto grossi. Ma come verificare che i risultati di inSPHector fossero corretti? La verifica della conservazione dell’energia che ci può dare un aiuto in tal senso, purtroppo non basta. Il debugging in certe situazioni non è la soluzione più adatta, vista la mole di dati che il programma riceve in input. Un approccio comodo è quello di rappresentare graficamente l’output del programma.


Come già detto inSPHector genera, durante la sua esecuzione, delle “fotografie” (chiamate snapshot) del sistema simulato. Rappresentando a video, uno dopo l’altro, tutti gli snapshot di una simulazione, si può verificare l’esattezza dei calcoli. Si analizzano gli spostamenti di una particella durante l’intero arco della simulazione, controllando che non sparisca o compia spostamenti troppo lunghi. Si è reso necessario uno strumento per la visualizzazione grafica degli snapshot. L’uso di applicazioni come GNUplot non è adatto, né sufficiente. StoB, ovvero “Snapshot to bitmap” risolve questi problemi. È uno strumento nato specificatamente per la generazione di immagini a partire da uno o più snapshot generati da inSPHector. Durante l’implementazione di StoB, abbiamo aggiunto delle caratteristiche che esulano dal puro scopo di controllare gli errori di inSPHector. Le immagini, infatti, ci possono dare importantissime informazioni qualitative su una simulazione, come la posizione degli agglomerati di particelle, la densità di questi agglomerati, etc. Ne descriveremo le caratteristiche più importanti, senza entrare nei dettagli.

4.2 Caratteristiche principali StoB è stato creato in Fortran 90, usando come modulo per l’accesso ai file, il modulo binary_m di cui abbiamo già parlato. Al programma serve necessariamente un file di testo (chiamato “stob.cfg”) che contiene informazioni sull’immagine da creare, come il colore di sottofondo, il tipo di immagine, etc.


Il programma ignora tutte le righe di questo file che sono vuote o che iniziano con il carattere “%”. L’ordine con cui i valori sono inseriti nel file è rilevante; per adesso infatti, StoB non esegue né un parsing né un syntax checking del file. Le immagini create hanno lo stesso nome dello snapshot che le genera, con in aggiunta il suffisso “.bmp”. Nel caso in cui venga generata in maniera automatica una serie di immagini, esse avranno prima del suffisso “.bmp” anche un numero progressivo di creazione. Le caratteristiche principali di StoB sono:

• facilità d’uso da linea di comando; • generazione di immagini in formato Windows bitmap; • determinazione automatica o manuale del dominio delle particelle

da rappresentare; • rappresentazione delle varie tipologie di particelle con colori e

simboli diversi; • supporto delle rotazioni e della prospettiva; • possibilità di creare delle immagini “quasi” fotografiche; • supporto di un linguaggio di script per le generazione di una serie di

immagini; • possibilità di specificare una palette dei colori.


4.3 Uso di StoB dalla linea di comando Digitando StoB, verrà stampato a video: Usage: [-i] <SnapshotFile> [-is|-s] <Script> [-l] <SnapshotList> [-ib] <BitmapFile> [-g|-h] -i Info on <SnapshotFile> -s Execute a script with name <Script> -is Total limits of snapshot in script with name <Script> -ib Print values used for <BitmapFile> -l Generate a image for all snapshots in file <SnapshotList> -g Generate a demo config file for stob -h Stob script 1.0 reference manual Stob version 1.0 prerelease 4, build 55 (30.01.01) Created by R.Favaretto, L.Moscardini, F.Sartoretto, G.Tormen University of Venice and University of Padova Per creare un’immagine basterà digitare: stob <SnapshotFile>. Se digitiamo : stob -i <SnapshotFile> otteniamo a video delle informazioni su <snapshotFile> come ad esempio: ----------------------------- INFORMATION ------------------------- Number of particles: 1472 gas : 1472 halo : 0 disk : 0 bulge: 0 stars: 0 Positions values min x: -0.9273847, max x: 0.9273847 min y: -0.9273847, max y: 0.9273847 min z: -0.9273847, max z: 0.9273847 Potential values min: -1.8926631, max: -0.9946011 EgySpec values min: 0.0500000, max: 0.0500000 Density values min: 0.0000000, max: 2.2248895 ------------------------------------------------------------------


Questo si rivela particolarmente utile per conoscere le dimensioni del dominio di <snapshotFile>, il numero di particelle e la loro tipologia. Con l’opzione -s possiamo decidere di eseguire uno script (verrà discusso più dettagliatamente nei prossimi paragrafi). Con -l <SnapshotList> possiamo far generare un’immagine per ogni snapshot elencato nel file di testo <SnapshotList>. L’opzione -g, StoB crea un test file Stob.cfg d’esempio. L’opzione –h, mostra a video una piccola guida all’uso del linguaggio di script incluso in Stob. È possibile inoltre, tramite l’opzione -ib <BitmapFile>, conoscere i parametri con cui è stata generata l’immagine <BitmapFile>.

4.4 Generazione di immagini in formato bmp Il programma genera immagini in formato bmp (Windows bitmap) a 16 milioni di colori a partire da snapshot generati da inSPHector. Questo formato grafico è stato scelto per:

• facilità di creazione di file di questo tipo; • supporto dei colori a 24 bit per pixel; • è universalmente riconosciuto dai programmi per la visualizzazione

grafica.


Uno svantaggio è la generazione di file grandi, poiché il formato bmp a 24 bit non supporta la compressione. In futuro verranno supportati altri tipi di file, ad esempio jpeg.

4.5 Determinazione automatica o manuale del dominio È possibile, tramite il file di configurazione, decidere se si vuole che venga determinato in maniera automatica il dominio della simulazione. In questo caso il programma adatta la posizione di tutte le particelle dello snapshot alla grandezza dell’immagine. Alternativamente è possibile impostare manualmente il dominio della simulazione. In questo caso si deve specificare il limite destro e sinistro lungo l’asse x, lungo l’asse y e lungo l’asse z del dominio. Questo è particolarmente utile se si vogliono effettuare degli zoom su uno snapshot. Nella figura 13, a sinistra vi è un’immagine creata con la determinazione automatica del dominio, mentre a destra con l’impostazione manuale (zoom del centro del cluster).


Figura 13: immagini generate da uno snapshot.

4.6 Rappresentazione delle particelle con colori e simboli diversi È possibile specificare nel file dei parametri il colore che ogni tipologia di particella deve avere nell’immagine. Ogni colore deve essere espresso tramite una terna RGB, ovvero da una terna di numeri compresi tra 0 e 255 che indicano rispettivamente il livello di rosso, il livello di verde ed il livello di blu presente nel colore da rappresentare. È possibile, inoltre, sfumare il colore di una particella in base a:

• distanza rispetto all’osservatore; • potenziale; • energia (solo per particelle SPH); • densità (solo per particelle SPH).

In questo caso si deve specificare il colore di partenza ed il colore d’arrivo.


Figura 14: rappresentazione di particelle con simboli diversi

È possibile rappresentare ogni tipo di particella con un simbolo diverso, in un elenco standard. Nella figura 14 sono state rappresentate le particelle di gas con un rombo, quelle di materia oscura con un cerchio.

4.7 Supporto delle rotazioni e della prospettiva Stob permette di ruotare le immagini rispetto agli assi cartesiani o rispetto al centro del dominio da rappresentare. Per fare questo si devono specificare, tramite una terna di valori, gli angoli di rotazione (in gradi) rispetto ai tre assi nel file di configurazione. Questa caratteristica è molto utile perché permette di poter vedere delle particelle che altrimenti sarebbero coperte da altre particelle.


Come ulteriore caratteristica, vi è la possibilità di creare delle immagini prospettiche.

4.8 Immagini “quasi” fotografiche Per creare, a partire da uno snapshot, delle immagini analoghe a quelle reali, abbiamo deciso di implementare un algoritmo che si basa sulla densità locale di ogni particella da rappresentare. In pratica per ogni particella si disegna un cerchio, di raggio inversamente proporzionale alla sua densità locale e colore di intensità direttamente proporzionale alla sua massa. Per fare questo abbiamo usato un oct-tree molto simile a quello usato per inSPHector.

Figura 15: esempio d’uso della modalità “foto”


Ogni cerchio può avere delle sfumature di tipo lineare, iperbolico, etc. Si può regolare l’immagine, variando il contrasto e la luminosità. Sono disponibili effetti fotografici per avere dei risultati congrui con diverse esigenze.

4.9 Il linguaggio di script

Abbiamo implementato in Stob un interprete di script, di formato proprietario, orientato alla creazione di sequenze di immagini. Facciamo un esempio: supponiamo che, dato lo snapshot con nome “snapshot” si vogliano creare 18 immagini partendo da una visione complessiva del dominio di simulazione, eseguendo uno zoom del 5 percento verso il centro e ruotando rispetto all’asse z di 20 gradi ad ogni immagine. Uno script che esegue queste operazioni è il seguente: repeat 18 create snapshot rotatez+ 20 zoomin 5 endrepeat Il linguaggio implementato non è case sensitive e consta di numerose funzioni: PRINT <string> Stampa a video una stringa


GLOBALROTATION Indica che le rotazioni devono essere eseguite rispetto il centro gli assi cartesiani

LOCALROTATION Indica che le rotazioni devono essere eseguite rispetto il centro dell’immagine SILENT <ON|OFF> Attiva o disattiva l’output a video SKIPOVERWRITE <ON|OFF> Se esiste già un’immagine con lo stesso

nome, con OFF viene sovrascritta, mentre con ON non viene creata

CREATE <snapshot> Crea un’immagine a partire dal file 'snapshot' REPEAT <intnum> Itera 'intnum' volte le istruzione che lo

seguono fino al corrispettivo ENDREPEAT ENDREPEAT Chiude un ciclo iniziato con REPEAT MINX <realnum> Imposta a 'realnum' il valore minimo per la coordinata x MAXX <realnum> Imposta a 'realnum' il valore massimo per la coordinata x MINY <realnum> Imposta a 'realnum' il valore minimo per la

coordinata y MAXY <realnum> Imposta a 'realnum' il valore massimo per la

coordinata y MINZ <realnum> Imposta a 'realnum' il valore minimo per la coordinata z


MAXZ <realnum> Imposta a 'realnum' il valore massimo per la

coordinata z MINX+ <realnum> Aggiunge 'realnum' al minimo valore per la

coordinata x MAXX+ <realnum> Aggiunge 'realnum' al massimo valore per la

coordinata x MINY+ <realnum> Aggiunge 'realnum' al minimo valore per la

coordinata y MAXY+ <realnum> Aggiunge 'realnum' al massimo valore per la

coordinata y MINZ+ <realnum> Aggiunge 'realnum' al minimo valore per la

coordinata z MAXZ+ <realnum> Aggiunge 'realnum' al massimo valore per la

coordinata z ROTATEX <realnum> Imposta a 'realnum' la rotazione rispetto l’asse x ROTATEY <realnum> Imposta a 'realnum' la rotazione rispetto l’asse y ROTATEZ <realnum> Imposta a 'realnum' la rotazione rispetto l’asse z ROTATEX+ <realnum> Aggiunge 'realnum' al corrente angolo di

rotazione lungo l’asse x ROTATEY+ <realnum> Aggiunge 'realnum' al corrente angolo di

rotazione lungo l’asse y


ROTATEZ+ <realnum> Aggiunge 'realnum' al corrente angolo di

rotazione lungo l’asse z MOVELEFT <realnum> Sposta il punto di vista a sinistra di 'realnum' MOVERIGHT <realnum> Sposta il punto di vista a destra di 'realnum' MOVEUP <realnum> Sposta il punto di vista verso l’alto di 'realnum' MOVEDOWN <realnum> Sposta il punto di vista verso il basso di 'realnum' MOVEFORWARDS <realnum> Sposta il punto di vista in avanti di 'realnum' MOVEBACKWARDS <realnum> Sposta il punto di vista indietro di 'realnum' ZOOMIN <realnum> Avvicinamento di 'realnum' percento del

punto di vista ZOOMOUT <realnum> Allontanamento di 'realnum' percento del

punto di vista

CAPITOLO 5. CONCLUSIONI 90

5. Conclusioni

Il lavoro presentato in questa tesi è stato sviluppato con l’intento di produrre un codice efficiente: sono state inserite una serie di ottimizzazioni che ci hanno permesso di migliorare notevolmente le prestazioni del programma rendendolo, nel caso in cui l’errore commesso nel calcolo dell’energia sia pressochè uguale, fino a otto volte più veloce rispetto a quello originale. Abbiamo implementato un BH-tree per ogni tipo di particella ed il calcolo dei vicini nel metodo SPH ora è effettuato usando un algoritmo che percorre il BH-tree stesso. Questa scelta ha migliorato l’efficienza rispetto all’algoritmo originario (algoritmo di Appeal) che si basa su un albero binario. È stato ottimizzato l’integratore temporale abbassando così il tempo necessario per effettuare l’integrazione dal 15% del tempo totale di una simulazione, a circa il 3%. È stato aggiunto il supporto delle coordinate comoventi e delle condizioni periodiche al contorno. Nel nuovo codice è possibile usare un criterio di apertura dell’oct-tree diverso da quello proposto da Barnes e Hut nel 1986. Questo nuovo criterio permette di avere un’accuratezza migliore nei calcoli laddove la distribuzione delle particelle è omogenea. Abbiamo reso il codice modulare, leggibile e ben strutturato, in modo da agevolare il più possibile le modifiche e la manutenzione. Sono state

CAPITOLO 5. CONCLUSIONI 91

ampiamente usate le “strutture” del Fortran 90. La memoria viene allocata dinamicamente. Abbiamo cercato di rendere l’uso del programma più semplice possibile inserendo la possibilità di leggere e scrivere i file sia in formato big endian che in formato little endian (a tale scopo è stata creato un modulo per l’accesso binario ai file). Il file dei parametri è stato completamente ristrutturato, prevedendo un controllo sintattico dello stesso in fase di avvio di una simulazione. Infine è stata creata una utility per la visualizzazione degli snapshot generati da inSPHector, che si è rilevata uno strumento utile per il debugging del programma, laddove era impensabile l’uso di strumenti alternativi come un debugger. Per quanto riguarda inSPHector, la fase successiva sarà quella di rendere il codice parallelo. L’utility stob va ottimizzata. Stiamo già valutando l’ipotesi di effettuare questi miglioramenti ottimizzando il codice attraverso l’uso delle librerie grafiche multi-piattaforma OpenGL. Un’interfaccia grafica, inoltre, ne renderebbe l’uso più facile, specialmente per la generazione di immagini di tipo fotografico. Questo lavoro è stato per me fonte di grande interesse e soddisfazione; probabilmente continuerà sotto forma di contratto con la facoltà di Astronomia dell’Università di Padova.

92

Ringraziamenti Desidero ringraziare il Prof. Flavio Sartoretto che ha dimostrato sempre interesse e disponibilità durante la realizzazione di questa tesi. Ringrazio i miei correlatori di Padova, Dott. Giuseppe Tormen e Dott. Lauro Moscardini, che mi hanno sopportato in questi mesi mentre cercavo di districarmi tra le difficili leggi che governano l’espansione dell’universo. Ringrazio i miei genitori, Rosa e Pio per la pazienza e per i finanziamenti che mi hanno elargito durante questi anni di università. Ringrazio la mia ragazza Francesca per le sere passate ad aspettarmi, mentre io restavo a lavorare alla realizzazione di questa tesi. Ringrazio tutti i miei amici in particolare Alessandro, Alessio, Filippo, Mauro e Stefano. Ringrazio tutti coloro che mi sono rimasti vicini e mi hanno supportato nel raggiungimento di questo obiettivo.

APPENDICE A. GLOSSARIO 93

A. Glossario

astrofisica: Disciplina che studia oggetti astronomici, come le stelle e le galassie.

barioni: Particelle pesanti che costituiscono la materia visibile.

Big Bang: La teoria del Big Bang postula un universo in espansione a partire da un punto infinitamente denso e caldo. L'istante iniziale è chiamato la "grande esplosione", il "Big Bang".

big endian: I byte di una variabile sono numerati in memoria da sinistra a destra, ovvero il byte più significativo ha indirizzo più basso rispetto al byte meno significativo (usato nelle CPU SPARC, Motorola e dei mainframe IBM).

carica: Numero quantico di una particella. Determina se la particella può partecipare a un dato processo di interazione. Una particella con carica elettrica ha interazioni elettriche, una particella con carica forte ha interazioni forti, etc.

carica elettrica: Il numero quantico che determina la partecipazione a interazioni elettromagnetiche.


CERN: Il principale laboratorio europeo internazionale, dotato di vari acceleratori, situato presso Ginevra, in Svizzera.

cluster (o ammasso): Famiglia di oggetti identici che in genere si suppone siano nati insieme da una struttura iniziale.

conservazione : Si dice che una quantità (carica elettrica, energia, quantità di moto) viene conservata quando è la stessa prima e dopo una reazione tra particelle.

conservazione della carica: Il fatto che la carica elettrica è conservata in ogni processo di trasformazione di un gruppo di particelle in un altro.

cosmologia: Lo studio della storia dell'universo.

costante di Hubble: L'universo osservabile è uniforme in tutte le direzioni, la sua velocità di espansione si può descrivere in termini di un solo parametro, la costante di Hubble, lentamente variabile nel tempo. La costante di Hubble è pari alla velocità media alla quale due regioni qualsiasi dell'universo si allontanano l'una dall'altra divisa per la distanza che le separa.

effetto Doppler: Effetto a causa del quale cambia la frequenza delle onde emesse da un oggetto in movimento. L'effetto permette di misurare la velocità di allontanamento delle galassie.

evento: Ciò che avviene quando due particelle collidono o una singola particella decade. A livello teorico si prevedono le probabilità di occorrenza di vari possibili eventi quando molte collisioni o decadimenti simili vengono studiati. Non è possibile prevedere il risultato di un singolo evento.

fotone : Il mediatore di forza delle interazioni elettromagnetiche.


gravitone : Il mediatore di forza delle interazioni gravitazionali; non è stato ancora osservato direttamente.

inflazione : Un periodo di espansione accelerata dell'universo, avvenuta pochi istanti dopo il Big Bang.

INFN: Istituto Nazionale di Fisica Nucleare, l'ente che in Italia coordina ed effettua ricerche sulla fisica delle particelle.

interazione : Processo in cui una particella decade o reagisce a una forza dovuta alla presenza di un'altra particella (come nelle collisioni).

interazione elettromagnetica: L'interazione dovuta alla carica elettrica; comprende le interazioni magnetiche.

interazione fondamentale: Nel Modello Standard le interazioni fondamentali sono le interazioni forte, elettromagnetica, debole e gravitazionale. La teoria prevede almeno un'ulteriore interazione fondamentale che sia responsabile delle masse delle particelle fondamentali. Cinque soli tipi di interazioni sono sufficienti e necessari per spiegare tutti i fenomeni fisici osservati.

interazione gravitazionale: L'interazione tra particelle dovuta alla loro massa/energia.

lensing gravitazionale: È il fenomeno di deflessione e distorsione della luce causato dalla gravità di un sistema, detto ”lente gravitazionale”, massiccio e con campo gravitazionale molto intenso (e.g. un ammasso di galassie).

little endian: I byte di una variabile sono numerati in memoria da destra a sinistra, ovvero il byte più significativo ha indirizzo più alto rispetto al byte meno significativo (usato nelle CPU Intel).


massa: Vedi "massa a riposo".

massa a riposo: La massa (m) in stato di quiete (o massa a riposo) di una particella è definita dall'energia della particella isolata (libera) in stato di quiete, divisa per il quadrato della velocità della luce.

materia oscura: Materia che esiste nello spazio ma non è visibile a noi perché non emette alcuna radiazione che la renda osservabile. Il movimento delle stelle attorno ai centri delle galassie cui appartengono presuppone che circa il 90% della materia di una galassia sia oscura. Sembra certo che ci sia materia oscura anche tra le galassie.

meccanica quantistica: Quadro concettuale che comprende le leggi fisiche che agiscono su scale piccolissime. Si basano sull'idea che la carica elettrica, la quantità di moto, il momento angolare, come anche le cariche, si danno in quantità discrete chiamate "quanti".

modulo: i moduli sono una caratteristica del Fortran 90 e sono molto simili alle unit del Pascal. Sono brani di codice che possono contenere funzioni e procedure, variabili globali, costanti e tipi derivati la cui visibilità agli altri moduli è limitata tramite l’uso delle parole chiave public e private.

nucleo: L'insieme di neutroni e protoni che forma il centro di un atomo.

parallasse: Effetto che risulta dalla differenza di aspetto di ogni oggetto visto da due punti differenti. Per analogia si indicano con questa parola delle quantità ben precise; la parallasse diurna (del Sole, ad esempio) è la differenza fra la direzione nella quale si vede contemporaneamente il Sole al tramonto in un certo luogo terrestre e quella nella quale si vede il Sole che sorge in un punto diametralmente opposto della Terra; è anche l'angolo sotto il quale si vede dal Sole il diametro terrestre. La parallasse annuale di una stella è l'angolo sotto il quale si vede da questa stella il semiasse maggiore dell'orbita della Terra intorno al Sole.


parsec: Unità di distanza, corrispondente alla distanza dal Sole di una stella, la cui parallasse annuale è di un secondo d'arco.

particella: Oggetto subatomico che ha massa e carica definite.

particella subatomica: Qualsiasi particella che sia più piccola di un atomo. private: parola chiave del Fortran 90 che permette di nascondere le funzioni, subroutine, variabili e costanti che la seguono, all’esterno del modulo in cui sono dichiarate. public: parola chiave del Fortran 90 che permette di rendere visibile le funzioni, subroutine, variabili e costanti che la seguono, anche all’esterno del modulo in cui sono dichiarate. quasar: Acronimo di ''Quasi Stellar Astronomical Radio Source”; si scrive anche QSS. Astri di massa paragonabile a quella delle galassie, molto concentrati, il cui spettro è caratterizzato da un eccesso di ultravioletto, un eccesso di infrarosso e righe di emissione del gas. redshift: Spostamento verso il rosso, per effetto Doppler, nella luce emessa da galassie lontane, dovuto alla loro recessione. snapshot: file, generati da inSPHector durante una simulazione, che rappresentano delle vere e proprie fotografie del sistema simulato. Permettono, inoltre, il riavvio di una simulazione. struttura: si veda type supernova: Evento catastrofico: una stella aumenta improvvisamente di luminosità, diventando visibile a grandissime distanze. Gli eventi di supernova sono usati per stimare la distanza delle galassie più lontane. Sono di due tipi: il tipo I, più spettacolare, riguarda stelle di popolazione anziana, il tipo II riguarda stelle di popolazione giovane.


type : parola chiave del Fortran 90 (simile alla struct del C e al record del Pascal) che permette di definire tipi derivati. La sintassi di base per definire un tipo derivato è la seguente: TYPE nomeNuovoTipo .

. [definizione variabile]

.

. END TYPE nomeNuovoTipo vento: Indica il movimento continuo di materia, che avviene a partire dal Sole, da una stella, da una galassia, diretto verso l'esterno (vento solare, vento stellare).

APPENDICE B. IL FILE DEI PARAMETRI 99

B. Il file dei parametri

Come già accennato nel paragrafo 3.4.3 il file dei parametri è un file di testo che viene letto dal programma all’inizio di ogni simulazione. In esso devono essere specificati tutti i parametri cosmologici necessari per definire il sistema fisico che si vuole evolvere e tutti i parametri che determinano le prestazioni del codice. Ogni parametro deve essere del tipo: nomeparametro <valore> Non è necessario che i parametri siano in un particolare ordine. Devono essere, però, tutti presenti. Il programma comunque effettua un syntax checking del file prima di far partire la simulazione e avverte l’utente se riscontra degli errori e/o delle incongruenze. Non è possibile, infatti, inserire parametri non previsti, come non è possibile inserire un parametro più di una volta. Il carattere di commento è il carattere “%”. Tutto ciò che lo segue, fino alla fine della riga, non viene preso in considerazione. Sono previste, inoltre, due istruzioni:

• %print <string> stampa a video <string> durante il parsing del file; • %sleep <sec> aspetta “sec” secondi durante il parsing del file.


Cercheremo qui di seguito di descrivere tutti i parametri che necessariamente devono essere inseriti nel file di configurazione. Parametro: InitialConditionsFile Esempio: InitialConditionsFile .\cluster.dat Descrizione: nome del file che contiene le condizioni iniziali per una

simulazione. Parametro: OutputFilesDir Esempio: OutputFilesDir .\output\ Descrizione: nome della directory di output del programma. In essa

verranno memorizzati tutti gli snapshot, i file statistici, i report, etc. generati durante una simulazione. Questa directory deve essere già presente, altrimenti il programma termina segnalando un errore.

Parametro: EnergyFile Esempio: EnergyFile energy.txt Descrizione: nome del file in cui verranno memorizzati, durante una

simulazione, i dati riguardanti l’energia del sistema. Parametro: InfoFile Esempio: InfoFile info.txt Descrizione: nome del file in cui verranno memorizzate, durante una

simulazione, informazioni riguardanti i passi temporali, gli snapshot generati, etc. È molto simile all’output a video.


Parametro: CpuFile Esempio: CpuFile cpu.txt Descrizione: nome del file in cui verranno memorizzate, durante una

simulazione, informazioni riguardanti il consumo di tempo cpu del programma.

Parametro: TimingsFile Esempio: TimingsFile timings.txt Descrizione: nome del file in cui verranno memorizzate, durante una

simulazione, varie statistiche sulle prestazioni del codice nel calcolo delle forze gravitazionali. Questo file viene creato solo se si è compilato il codice con la direttiva di precompilazione DIAG.

Parametro: DefaultSaveFormat Esempio: DefaultSaveFormat 0 Descrizione: formato con cui devono essere memorizzati i file binari

generati durante una simulazione. Questo parametro può assumere tre valori: • 0: i file vengono memorizzati nel formato tipico della macchina in cui sta girando il programma. • 1: il file viene memorizzato in formato little endian • 2: il file viene memorizzato in formato big endian

Parametro: SnapshotPrefix Esempio: SnapshotPrefix snapshot_ Descrizione: nome di ogni snapshot è composto da questa stringa con in

aggiunta un numero progressivo a tre cifre.


Parametro: RestartFile Esempio: RestartFile restart Descrizione: nome del file di riavvio di una simulazione. Se non si

specifica alcun valore, esso non viene mai creato, altrimenti viene creato sia ad ogni intervallo di tempo (si veda il parametro TimeBetRestartFile), sia alla fine della simulazione.

Parametro: SnapshotTimeList Esempio: SnapshotTimeList list Descrizione: nome del file che contiene la lista dei tempi in cui generare

gli snapshot. Se non viene specificato niente vengono creati gli snapshot in base al parametro NumberOfSnapshot. I valori elencati in questo file non occorre che siano ordinati per tempo.

Parametro: MaxAllowedCPU Esempio: MaxAllowedCPU 1000 Descrizione: esprime, in secondi, il tempo massimo d’esecuzione del

programma. Se si è prevista la scrittura di un file di riavvio (si veda il parametro restartFile), InSPHector superato l’85% di questo valore inizia la procedura di arresto della simulazione. Il restante 15% serve per avere la sicurezza che venga memorizzato il file di riavvio e l’eventuale snapshot finale. Nel caso in cui non sia prevista la scrittura del file di riavvio, il programma superato il 95% del tempo specificato inizierà la procedura di arresto. Il restante 5% servirà per la scrittura dell’eventuale snapshot finale. Se il valore specificato è zero, il programma non ha limiti di tempo per la sua esecuzione.


Parametro: ResubmitScript Esempio: ResubmitScript script Descrizione: file che può essere eseguito nel caso in cui la simulazione

venga interrotta poiché il programma ha superato il tempo cpu a sua disposizione (si veda il parametro MaxAllowedCPU).

Parametro: InitialConditionFileFormat Esempio: InitialConditionFileFormat 3 Descrizione: tipologia del file delle condizioni iniziali che il programma

deve leggere per avviare una simulazione. Parametro: TimestepType Esempio: TimestepType 1 Descrizione: può assumere valore 0 o 1. Questo parametro seleziona il

criterio per il calcolo dei passi temporali delle particelle. Se vale 1 ogni passo temporale è proporzionale a a/1 dove a è l’accelerazione della particella presa in esame. Se vale 0 è proporzionale a a/1 .

Parametro: StartingTime Esempio: StartingTime 2.0 Descrizione: tempo di inizio della simulazione.


Parametro: MaximumTime Esempio: MaximumTime 4.0 Descrizione: tempo finale della simulazione. Se

MaximumTime<StartingTime allora essi vengono interpretati come valori di redshift. Di conseguenza viene abilitato l’uso delle coordinate comoventi.

Parametro: Omega0 Esempio: Omega0 0.3 Descrizione: parametro cosmologico della densità critica a redshift uguale

a zero. Rilevante solo nel caso di coordinate comoventi. Parametro: OmegaLambda Esempio: OmegaLambda 0 Descrizione: parametro cosmologico della densità critica nel vuoto a

redshift uguale a zero. Rilevante solo nel caso di coordinate comoventi. Per un universo piatto Omega0+OmegaLambda=1.

Parametro: OmegaBaryon Esempio: OmegaBaryon 0 Descrizione: parametro cosmologico della densità critica della materia

barionica a redshift uguale a zero. Rilevante solo nel caso di coordinate comoventi.


Parametro: HubbleConstant Esempio: HubbleConstant 1.0 Descrizione: valore della costante di Hubble espressa in unità di

11100 −−⋅ MpcsKm . Necessaria nelle simulazioni in cui sono necessarie conversioni al sistema cgs.

Parametro: BoxSize Esempio: BoxSize 0 Descrizione: grandezza del cubo contenente tutte le particelle del sistema

da simulare. Rilevante solo se il codice è stato compilato con l’opzione di precompilazione PERIODIC. Le coordinate delle particelle devono essere comprese nell’intervallo [0,BoxSize].

Parametro: NumberOfSnapshot Esempio: NumberOfSnapshot 10 Descrizione: numero di snapshot che si desidera vengano scritti durante una simulazione. Parametro: TimeOfFirstSnapshot Esempio: TimeOfFirstSnapshot 0 Descrizione: tempo dopo il quale iniziare a scrivere gli snapshot.


Parametro: TimeBetRestartFile Esempio: TimeBetRestartFile 3600.0 Descrizione: intervallo di tempo, espresso in secondi, passato il quale il

programma deve memorizzare un file di riavvio (vedere RestartFile) . Questo parametro rappresenta una misura precauzionale per i limitare i danni dovuti a problemi hardware o software del sistema durante una simulazione. Se vale zero, nessun file di riavvio, tranne l’ultimo, viene creato.

Parametro: TimeBetStatistics Esempio: TimeBetStatistics 0.05 Descrizione: intervallo di tempo, espresso in secondi, passato il quale il

programma deve memorizzare delle informazioni, nel file specificato al parametro EnergyFile, riguardanti l’energia del sistema.

Parametro: DeltatTollerance1 Esempio: DeltatTollerance1 0.2 Descrizione: nel caso in cui TimestepType sia uguale a 0 nel codice

abbiamo usato un criterio alternativo per il calcolo dei passi

temporali. In questo caso a/2ηε=∆t dove

erance1DeltatToll=η e ε è la lunghezza del softening gravitazionale.


Parametro: DeltatTollerance2 Esempio: DeltatTollerance2 0.05 Descrizione: usato solo se TimestepType è uguale a 1. Rappresenta una

frazione della velocità tipica delle particelle ed è usato per il calcolo dei passi temporali . In questo caso

a/erance2DeltatTollt =∆ .

Parametro: CourantFactor Esempio: CourantFactor 0.1 Descrizione: valore di Courant per il calcolo dei passi temporali delle particelle SPH Parametro: MaxTimestepLength Esempio: MaxTimestepLength 0 Descrizione: massimo passo temporale necessario allo spostamento di una particella. Parametro: MinTimestepLength Esempio: MinTimestepLength 0 Descrizione: minimo passo temporale necessario allo spostamento di una particella.


Parametro: TreeOpeningCriterionType Esempio: TreeOpeningCriterionType 1 Descrizione: criterio di apertura usato per la visita del BH-albero. Se

zero, viene usato il criterio standard di Barnes-Hut, mentre se vale uno, viene usato il criterio descritto da V.Springel, N.Yoshida e S.D.M.White (2000).

Parametro: TreeOpeningAngle Esempio: TreeOpeningAngle 0.8 Descrizione: parametro di accuratezza per l’algoritmo di barnes-Hut. Se

TreeOpeningCriterionType vale 1, questo valore viene usato solo per il primo calcolo delle forze gravitazionali. Successivamente si userà il valore assegnato a TreeForceAcc.

Parametro: TreeForceAcc Esempio: TreeForceAcc 0.02 Descrizione: parametro di acuratezza per il criterio di apertura dell’albero

introdotto da V.Springel, N.Yoshida e S.D.M.White (2000). Rilevante solo se il parametro TreeOpeningCriterionType vale uno.

Parametro: TreeUpdatingLimit Esempio: TreeUpdatingLimit 0.1 Descrizione: questo parametro indica quanti calcoli di forza rispetto il

numero di particelle del sistema si devono fare prima di ricostruire tutto il BH-albero. Ad esempio se N è il numero di particelle del sistema e TreeUpdatingLimit vale N/10, allora solamente dopo N*(N/10) calcoli di forze l’albero


viene interamente ricostruito. Se questo parametro vale 0, l’albero viene ricostruito ad ogni passo temporale.

Parametro: MovingNodesLimit Esempio: MovingNodesLimit 0.05 Descrizione: questo parametro protegge l’algoritmo da rapidi

cambiamenti di parti dell’albero anche nel caso in cui esso non sia aggiornato di frequente. Sia i un generico nodo: se il suo centro di massa si è mosso più di un fattore MovingNodesLimit volte la lunghezza del suo lato dall’ultima ricostruzione dell’albero, allora i viene completamente aggiornato, ossia i suoi momenti multipolo ed il suo centro di massa sono ricalcolati.

Parametro: SPHNeighbours Esempio: SPHNeighbours 20 Descrizione: numero di vicini che si desiderano per ogni particella SPH. Parametro: ArtificialViscosity Esempio: ArtificialViscosity 0.75 Descrizione: parametro di viscosità artificiale. Parametro: StartingGasTemperature Esempio: StartingGasTemperature 0 Descrizione: temperatura iniziale, espressa in Kelvin, del gas. La

temperatura delle particelle di gas assume questo valore solo se esso è maggiore di zero e se la loro temperatura vale zero.


Parametro: MinimumGasTemperature Esempio: MinimumGasTemperature 0 Descrizione: temperatura minima che possono avere la particelle di gas. Parametro: cgsLengthUnit Esempio: cgsLengthUnit 1.0 Descrizione: assegna le unità di misura del programma in cm/h dove

110 100 −−⋅⋅= MpcsKmhH

Parametro: cgsMassUnit Esempio: cgsMassUnit 1.0 Descrizione: assegna le unità di misura del programma in g/h dove

110 100 −−⋅⋅= MpcsKmhH

Parametro: cgsVelocityhUnit Esempio: cgsVelocityhUnit 1.0 Descrizione: assegna le unità di misura del programma in cm/s. Parametro: CodeGravitationalConstant Esempio: CodeGravitationalConstant 1.0 Descrizione: valore della costante gravitazionale riferita alle unità di

misura scelte. Se vale 0 il programma la calcola automaticamente.


Parametro: MinSPHtoGravSoftRatio Esempio: MinSPHtoGravSoftRatio 0.0 Descrizione: minima lunghezza di smussamento per le particelle di gas in

termini di lunghezza di smussamento gravitazionale. Parametro: GasSoftening Esempio: GasSoftening 0.04 Descrizione: softening per le particelle di gas in unità interne al

programma. Parametro: HaloSoftening Esempio: HaloSoftening 0 Descrizione: softening per le particelle di materia oscura in unità interne

al programma. Parametro: DiskSoftening Esempio: DiskSoftening 0 Descrizione: softening per i “dischi stellari” in unità interne al

programma. Parametro: BulgeSoftening Esempio: BulgeSoftening 0 Descrizione: softening per i “nuclei stellari” (bulge) in unità interne al

programma.


Esempio: StarsSoftening 0 Descrizione: softening per le stelle in unità interne al programma. Parametro: MaxPhysGasSoftening Esempio: MaxPhysGasSoftening 0.04 Descrizione: massimo softening fisico per le particelle di gas. Parametro: MaxPhysHaloSoftening Esempio: MaxPhysHaloSoftening 0 Descrizione: massimo softening fisico per le particelle di materia oscura. Parametro: MaxPhysDiskSoftening Esempio: MaxPhysDiskSoftening 0 Descrizione: massimo softening fisico per i “dischi stellari”. Parametro: MaxPhysBulgeSoftening Esempio: MaxPhysBulgeSoftening 0 Descrizione: massimo softening fisico per i “bulge stellari”. Parametro: MaxPhysStarsSoftening Esempio: MaxPhysStarsSoftening 0 Descrizione: massimo softening fisico per le stelle.

APPENDICE C. ALBERI ROSSO-NERI 113

C. Alberi rosso-neri Nel nostro codice abbiamo usato un albero rosso-nero per memorizzare gli indice delle particelle, ordinate rispetto al tempo di spostamento. La sua implementazione è stata fatta attraverso l’ausilio di un array di dimensione pari al numero di particelle, N, presenti nel sistema. Questo non ci sembra limitativo, dato che l’albero, durante una simulazione, ha sempre N nodi. Un albero rosso-nero è un albero binario di ricerca. Per ogni nodo un campo binario contiene il suo colore che può essere rosso o nero. Sia N il numero di nodi dell’albero. Attraverso precise regole di colorazione dei nodi, su qualsiasi cammino dalla radice ad una foglia, nessun cammino di un albero rosso-nero risulta lungo più del doppio di qualsiasi altro. Diremo quindi che l’albero è approssimativamente bilanciato. Una conseguenza diretta di questo, consiste nel fatto che gli algoritmi di ricerca e di inserimento hanno costo )(log NO nel caso pessimo.


Ricordiamo che in un albero binario di ricerca, la ricerca e l’inserimento di un nodo costano O(N) nel caso in cui l’albero sia sbilanciato (caso pessimo). Diremo che un albero binario di ricerca è un albero rosso-nero se soddisfa le seguenti proprietà:

• ogni nodo può essere solamente rosso o nero; • le foglie sono tutte nere; • se un nodo è rosso, entrambi i figli sono neri; • ogni cammino a partire da un nodo ad una qualsiasi delle foglie

contiene lo stesso numeri di nodi neri.

20

412918NIL

3515

NILNIL NIL NILNIL NIL

Figura 16: esempio di albero rosso-nero

Dopo aver inserito un nuovo elemento all’interno dell’albero (visto come albero binario di ricerca) bisogna procedere con una serie di operazioni in modo tale da riottenere un albero rosso-nero.


Descriveremo tutte queste operazioni senza mostrarne l’implementazione; per approfondire l’argomento si può consultare ad esempio [Cormen, Leiserson e Rivest, 1994]. Rotazione Le operazioni di inserimento e cancellazione di un nodo su un albero rosso-nero con N nodi impiegano tempo )(log NO . Tuttavia, dato che esse modificano la struttura dell’albero stesso, è possibile che l’albero prodotto violi le proprietà degli alberi rosso-neri. Per ripristinare queste proprietà è necessario cambiare il colore a qualche nodo dell’albero ed eventualmente modificarne la struttura. Le operazioni preposte a questo scopo prendono il nome di rotazioni. Esse possono essere effettuate verso sinistra oppure verso destra. Quando effettuiamo una rotazione a sinistra su un nodo x, assumiamo che il suo figlio destro y sia diverso da NIL. La rotazione a sinistra pone y come nuova radice del sottoalbero, x come figlio sinistro di y e il figlio di y come figlio destro di x. La rotazione a destra è simile alla rotazione a sinistra. Entrambe hanno costo O(1). Inserimento L’inserimento di un nodo in un albero rosso-nero con N nodi può essere realizzato in tempo )(log NO . Inizialmente si procede nell’inserimento del nuovo nodo x colorandolo di rosso e considerando l’albero rosso-nero come un albero binario di ricerca. Successivamente, al fine di garantire le proprietà degli alberi rosso-neri, si procede ad una serie di ricolorazioni e rotazioni dei nodi.


Ricerca Tipica degli alberi binari di ricerca. Nel caso peggiore richiede tempo

)(log NO poiché gli alberi rosso-neri sono approssimativamente bilanciati. Cancellazione Come per le altre operazioni di base su un albero rosso-nero con N nodi, la cancellazione di un nodo richiede tempo )(log NO . Poiché essa cambia la struttura dell’albero, è necessario, attraverso l’ausilio delle rotazioni, ripristinare le proprietà “rosso-nere” dell’albero.

APPENDICE D. IL MODULO BINARY_M 117

D. Il modulo binary_m Come già accennato, inSPHector è in grado di riconoscere automaticamente se i file delle condizioni iniziali, sono in formato big endian o in formato little endian. Questa caratteristica è stata inserita nel codice, poiché a volte risulta necessario terminare una simulazione su di un computer con architettura diversa da quella del computer in cui è stata iniziata. Poichè il Fortran 90 standard non permette di poter gestire i file binari in maniera sequenziale, si è resa necessaria la creazione di un modulo che permettesse di fare questo. A livello implementativo, abbiamo usato il supporto che fornisce il Fortran ai file sequenziali, “bufferizzando” la lettura e la scrittura su file. Quello che abbiamo ottenuto è stato un set completo di routine che assomigliano alle istruzioni standard del Fortran (come Open, Close, Read e Write) e che quindi sono molto facili da usare, agevolando moltissimo l’utente nello scrivere codice “orientato” ai file binari. In aggiunta abbiamo inserito delle funzioni di utilità generale che descriveremo meglio più avanti in questa appendice.


Come già accennato, l'uso delle routine di questo modulo è molto simile all'uso delle istruzioni standard del Fortran 90; le differenze sono:

• ogni subroutine ha prefisso bin; • ogni parametro ha nome ottenuto dallo standard Fortran 90,

aggiungendo il prefisso bin; • si deve aggiungere call per la chiamata della subroutine.

! Apre il file in scrittura

call binOpen(fd,binfile="strings.bin",binstatus="replace", &

binaction="write")

call binWrite(fd,10) ! scrive 10 nel file “strings.bin”

call binClose(fd) ! chiude il file

! Apre il file in lettura

call binOpen(fd,binfile="strings.bin",binstatus="old", &

binaction="read")

call binRead(fd,intero) ! legge un intero

call binClose(fd) ! chiude il file

print *,intero ! stampa 10 a video

Esempio d’uso di alcune routine del modulo

L'elenco delle routine richiamabili è il seguente:

• binOpen; • binClose; • binFlush; • binRewind; • binConversion; • binWrite; • binRead.


Esse verranno descritte più approfonditamente qui di seguito. Subroutine: binOpen Descrizione: apre un file binario Parametri: BINUNIT (obbligatorio): numero intero, compreso tra 90 e 99, da associare

al file; BINFILE (obbligatorio): nome del file da aprire; BINACTION (facoltativo): può essere scelto tra “READ”, “WRITE” e

“READWRITE” e rappresenta la modalità con cui si intende aprire il file: lettura, scrittura e lettura-scrittura;

BINSTATUS (facoltativo): descrive la natura del file da aprire. Può essere uno dei seguenti:

• “OLD” : apre un file esistente. Se non esiste viene prodotto un errore;

• “NEW”: il file da aprire non deve esistere. Se esiste viene prodotto un errore;

• “SCRATCH”: il file viene aperto e una volta chiuso, viene automaticamente cancellato;

• “REPLACE”: il file aperto rimpiazza il file (se esiste) con lo stesso nome;

• “UNKNOWN” (Default). Il sistema a run-time dapprima cerca di aprire il file con binstatus uguale a “OLD”, e poi, se non esiste, lo apre con “NEW”.

BINIOSTAT (facoltativo): intero che vale 0 se non ci sono stati errori; BINPOSITION (facoltativo): può essere scelto tra “REWIND” e

“APPEND“. Con il primo il puntatore al file è all’inizio del file stesso, mentre con il secondo è alla fine.


Subroutine: binClose Descrizione: chiude un file binario Parametri: BINUNIT (obbligatorio): numero intero che rappresenta il file da chiudere; BINIOSTAT (facoltativo): intero che vale 0 se non ci sono stati errori; BINSTATUS (facoltativo): può essere “KEEP” (default) oppure

“DELETE”. Nel secondo caso il file viene eliminato dopo la chiusura.

Subroutine: binFlush Descrizione: svuota il buffer Parametri: BINUNIT (obbligatorio): numero intero che rappresenta il file il cui buffer

deve essere svuotato; BINIOSTAT (facoltativo): intero che vale 0 se non ci sono stati errori. Subroutine: binRewind Descrizione: ritorna all’inizio del file Parametri: BINUNIT (obbligatorio): numero intero che rappresenta il file. Subroutine: binConversion Descrizione: svuota il buffer Parametri t: BINUNIT (obbligatorio): numero intero che rappresenta il file; BINVAL (obbligatorio): se .true. ed il file aperto è in formato big endian e

il sistema lavora in rappresentazione little endian, i dati letti o scritti vengono convertiti in big endian. Se .true. ed il file aperto è in formato


little endian e il sistema lavora in rappresentazione big endian, i dati letti vengono convertiti in little endian.

Subroutine: binWrite Descrizione: scrive su disco uno o più byte Parametri: BINUNIT (obbligatorio): numero intero che rappresenta il file; BINVAL: (obbligatorio): variabile il cui contenuto è da scrivere su file.

Essa può essere di tipo: character(*), integer(1), integer(2), integer(4), real(4), real(8) e logical(4);

BINIOSTAT (facoltativo): intero che vale 0 se non ci sono stati errori. oppure: BINUNIT (obbligatorio): numero intero che rappresenta il file; BINARRAY (obbligatorio): array il cui contenuto va scritto su file. Ogni

elemento dell’array può essere di tipo: character(*), integer(1), integer(2), integer(4), real(4), real(8) e logical(4);

BINELEMENTS (obbligatorio): numero di elementi da scrivere; BINIOSTAT (facoltativo): intero che vale 0 se non ci sono stati errori. Subroutine: binRead Descrizione: legge da disco uno o più byte Parametri: analogo a binWrite.


Nel modulo sono incluse delle funzioni aggiuntive e di utilità: Function: machineArch Descrizione: determina l’architettura della macchina su cui sta girando il programma Parametri: nessuno. Valori di ritorno: 0 se il sistema è little endian; 1 se il sistema è big endian; -1 se il sistema ha una rappresentazione dei dati sconosciuta. Function: typeOfNotFormattedFile Descrizione: determina la tipologia di un file non formattato Parametri: fileName (obbligatorio): nome del file. Valori di ritorno: 0 se il file è in formato little endian; 1 se il file è in formato big endian; -1 se il file è di formato sconosciuto oppure non è un file valido; -2 se la funzione non riesce a capire se il file è big o little endian; -3 se il file non è stato trovato. Function: fileLength Descrizione: restituisce la lunghezza in byte di un file Parametri: fileName (obbligatorio): nome del file. Valori di ritorno: lunghezza (in byte) del file; -1 in caso di errore.

BIBLIOGRAFIA 123

Bibliografia [1] F. SARTORETTO, M. PUTTI, Fortran per applicazioni numeriche, Edizioni Libreria Progetto Padova, 2000. [2] J. F. NAVARRO, S. D. M. WHITE, Simulations of dissipative galaxy formation in hierarchically clustering uiversis - I. Tests of the code, Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, Physics Departement, University of Durham, Durham DH1 3LE, 1993. [3] T. J. MARTIN, F. R. PEARCE, P. A. THOMAS, An Owner’s Guide to Smoothed Particle Hydrodynamics, Astronomy Centre, Sussex University, Falmer, Brighton, BN1 9 QH, 1994. [4] G. TORMEN, Origine ed evoluzione dell’universo, Dipartimento di Astronomia dell’Università di Padova, 2000. [5] E. BERTSCHINGER, Cosmic Structure Formation, Institute for advanced Study, Princeton, NJ 08540 USA, Departement of Physics, MIT 6-207, Cambridge, MA 02139, USA, 1993. [6] F. A. RASIO, Particle Methods in Astrophysical Fluid Dynamics, Department of Physics, MIT, Cambridge, MA 02139, USA, 1999. [7] V. SPRINGEL, N. YOSHIDA, S. D. M. WHITE, GADGET: A code for collisionless and gasdynamical cosmological simulations, Max-Planck-Insitut fur Astrophysik, Karl-Schwarzschild- Strabe 1, 85740 Garching bei

BIBLIOGRAFIA 124

Munchen, Germany, Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA, 1993. [8] A. GULLI, Formulazioni parallele dei metodi di simulazione N-Body, tesi di laurea, Università degli Studi di Pisa, 1997. [9] T. H. CORMEN, C. E. LEISERSON, R. L. RIVEST, Introduzione agli algoritmi, Jackson libri, Milano, vol. 1, 1994. [10] J. D. BARROW, Le origini dell’universo, Sansoni Editore, Milano, 1995. [11] D. W. SCIAMA, Cosmologia Moderna, Arnoldo Mondadori Editore Spa, Milano, 1999. [12] S. RIMMER, Bit Mapped Graphics, McGraw-Hill Libri Italia srl, Milano, 1993. [13] J. E. BARNES, Treecode Guide, Institute for Astronomy, University of Hawai’i, 2680 Woodlawn Drive, Honolulu, Hawai’i 96822. [14] J. E. BARNES, P. HUT, The use of Supercomputer in Stellar Dynamics, eds. P. Hut e S. Mc Millan, Berlin: Springer-Verlag, 1986. [15] W. BENZ, Smoothed particle hydrodynamics: A review, Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA, 1993. [16] S. J. CHAPMAN, Fortran 90/95 for Scientists and Engineers, McGraw-Hill, Boston, 1996. [17] C. C. DYER, S. S. PETER, Softening in N-body simulations of collisionless system, Astrophysical Jurnal Supplement, 409, 60-67, 1993. [18] L. HERNQUIST, F. R. BOUCHET, Y. SUTO, Application of the Ewald method to cosmological N-body simulations, Astrophysical Journal Supplement, vol. 75, 231-240, 1991.

BIBLIOGRAFIA 125

[19] L. HERNQUIST, N. KATZ, TREESPH - A unification of SPH with the hierarchical tree method, Astrophysical Journal Supplement, vol. 70, n.2, 419-446, 1989. [20] J. J. MONAGHAN, Particle Methods for Hydrodynamics, Computer Physics Reports 3, 71-124, 1985.

INDICE DELLE FIGURE 126

Indice delle figure Figura 1: condizioni iniziali della simulazione .............................................5

Figura 2: risultato di una simulazione idrodinamica .....................................7

Figura 3: funzione kernel usata nel codice ..................................................13

Figura 4: esempio di oct-tree .......................................................................19

Figura 5: logica del programma ..................................................................23

Figura 6: esempio di multi-albero e relativa rappresentazione con array ...39

Figura 7: approssimazione di Burnes-Hut ...................................................43

Figura 8: schema dei moduli e delle relative inclusioni ..............................55

Figura 9: Simulazione 1: tempi di calcolo...................................................72

Figura 10: Simulazione 1: conservazione dell’energia. ..............................74

Figura 11: Simulazione 2: tempi di calcolo.................................................75

Figura 12: Simulazione 2: conservazione dell’energia ...............................76

Figura 13: immagini generate da uno snapshot. ..........................................83

Figura 14: rappresentazione di particelle con simboli diversi.....................84

Figura 15: esempio d’uso della modalità “foto”..........................................85

Figura 16: esempio di albero rosso-nero ...................................................114

INDICE ANALITICO

127

Indice analitico alberi rosso-neri.................................................................................115; 116 albero ...... 5; 18-20; 36-42; 44; 50-<54; 64; 90; 108; 109; 113; 114; 115; 116 arresto ..................................................................................................30; 102 arresto forzato ..............................................................................................30 barioni..........................................................................................................93 Benz.............................................................................................................21 big endian ..............................................6; 30; 32; 91; 93; 101; 117; 120; 122 cancellazione ...............................................................................43; 115; 116 caratteristica.........................................................5; 8; 10; 46; 84; 85; 96; 117 cella ............................................................................................19; 44; 51; 52 central leap-frog...........................................................................................21 centro di massa ......................................................18; 19; 20; 35; 40; 42; 109 condizioni periodiche ............................................................6; 22; 46; 57; 90 coordinate comoventi ....................................................6; 22; 33; 45; 90; 104 coordinate fisiche .........................................................................................33 cosmologia...........................................................................................4; 8; 94 densità ..................5; 7; 11; 12; 14; 18; 20; 35; 43; 45; 53; 61; 78; 83; 85; 104 dominio ....................................................................18; 46; 79; 81; 82; 84; 86 efficienza ...............................................................................6; 71; 75; 76; 90 file dei parametri..............................16; 23; 25; 28; 30; 32; 33; 41; 83; 91; 99 file di configurazione .......................................26; 28; 48; 53; 56; 82; 84; 100 Fortran 90 ........................6; 21; 24; 25; 34; 78; 91; 96; 97; 98; 117; 118; 124 forze ...........4; 5; 8; 9; 10; 14; 17; 22; 40; 41; 42; 43; 44; 48; 56; 61; 101; 108

INDICE ANALITICO

128

gravitazione .....................................................................................16; 17; 45 immagine .................................................................78; 80; 81; 82; 83; 86; 87 integratore temporale ...................................................................6; 21; 47; 90 kernel.............................................................................10; 12; 13; 59; 64; 67 leggi idrodinamiche ...............................................................................21; 35 linea di comando........................................................................25; 26; 79; 80 linguaggio ..................................................................................24; 79; 81; 86 little endian..................................................6; 32; 91; 95; 101; 117; 120; 122 lunghezza di smussamento ..........................................15; 35; 48; 52; 53; 111 modifiche strutturali ..............................................................................22; 24 modulo ................................. 24; 25; 32; 56-63; 78; 91; 96; 97; 117; 118; 122 oct-tree .......................6; 18; 19; 20; 25; 29; 36; 37; 41; 42; 51; 57; 75; 85; 90 palette ..........................................................................................................79 particella 4; 6; 10-25; 34-38; 42-54; 57; 61; 78; 83-85; 90-97; 103; 107; 109 particelle barioniche ..........................................................5; 8; 35; 38; 53; 62 particelle di materia oscura..................................................5; 8; 38; 111; 112 passo temporale ...................................................................38; 103; 107; 109 pressione ......................................................................................5; 13; 14; 35 print ..........................................................................................24; 29; 99; 118 riavvio ........................................................................24; 26; 30; 97; 102; 106 script ..................................................................................79; 80; 81; 86; 103 simulazione ....... 4; 5; 7; 18-61; 71-78; 82; 86; 90; 91; 97-106; 113; 117; 124 sleep .......................................................................................................29; 99 Smoothed Particle Hydrodynamics .......................................................9; 123 snapshot ..... 24; 26; 27; 30; 31; 33; 60; 78-86; 87; 91; 97; 100; 101; 102; 105 spline ............................................................................................................12 vicini ....................................6; 16; 35; 51; 52; 53; 62; 72; 73; 76; 90; 92; 109

Dipartimento di Fisica e Astronomia - Indiceinsieme di equazioni differenziali ordinarie e le forze...

Documents

Transcript of Dipartimento di Fisica e Astronomia - Indiceinsieme di equazioni differenziali ordinarie e le forze...