Applicazioni di Calcolo
Parallelo in Astrofisica
Claudio Gheller
CINECA
[email protected]
Firenze, 10-11 Giugno 2003, C. Gheller
Astronomia e Computer: un legame molto stretto.
DALLE
OSSERVAZIONI
ALLA TEORIA
L’uso del calcolatore è indispensabile in ogni passo della ricerca
astronomica:
• Sviluppo e progettazione di nuovi strumenti di osservazione
• Controllo e gestione delle strumentazioni (telescopi, satelliti)
• Acquisizione dei dati, telemetria
• Data management e compressione
• Riduzione e analisi dei dati
• Elaborazione dei risultati
• Realizzazione di modelli teorici
• Confronto tra osservazioni e modelli teorici
Firenze, 10-11 Giugno 2003, C. Gheller
Simulazioni Numeriche
Le simulazioni numeriche, sono la
realizzazione dei modelli teorici
• Descrizione del modello
• Condizioni iniziali
• Equazioni evolutive
• Rappresentazione dello spazio
• Approssimazione del modello
• Algoritmi di soluzione delle equazioni
• Ricostruzione degli osservabili
• Confronto con le osservazioni
Firenze, 10-11 Giugno 2003, C. Gheller
Simulazioni numeriche in astrofisica
I campi applicativi sono
innumerevoli:
Il sistema solare
Il sole, nascita ed evoluzione dei
pianeti…
Le stelle
Struttura, evoluzione, ammassi…
Le galassie
Nascita, evoluzione…
La cosmologia
Proprietà ed evoluzione dell’Universo
Firenze, 10-11 Giugno 2003, C. Gheller
Cosmological Simulations
Cosmological simulations are
challenging applications since:
There are no special
symmetries.
Full 3-D simulations (periodic
boundary conditions);
No special shapes of the forming
structures.
Two different kind of matter to
be treated (different models –
different alghoritms).
Dark Matter
Baryons
Must describe a
representative fraction of
the universe.
10 – 100 Mpc;
i.e. 1020 – 1021 km.
Cosmological structures are
influenced by phenomena on a
n enormous range of scales
Star forming regions (1 pc);
Galaxies (10 kpc);
Clusters of galaxies (1 Mpc);
Pressure forces (10 Mpc);
Gravity (100 Mpc).
DYNAMICAL RANGE
=
109 !!!
N-Body;
Hydrodynamics.
What is a Cosmological
Simulation ???
A cosmological simulation
has the object of describing
the birth and the evolution of
the cosmological structures
(galaxies, galaxy clusters,
filaments)
following
the
dynamics
of
all
the
components:
baryons
(typically fluid description),
dark matter and stars (Nbody)…
Firenze, 10-11 Giugno 2003, C. Gheller
Applicazioni Cosmologiche
Viene simulato un volume di dimensioni
confrontabili all’intero universo (1/100)
Cosa possiamo imparare ???
Fisica fondamentale e universo primordiale
Interazioni fondamentali, Big Bang, Inflazione,
transizioni di fase…
COBE mission:
map of the infrared sky
Dipole map
Dipole subtracted map
La più recente immagine
del Big Bang (MAP)
Milky way subtracted map
T=2.73 K, dT/T = 10-5
Resolution 10°
Firenze, 10-11 Giugno 2003, C. Gheller
Cosa possiamo imparare???
Struttura e proprietà delle galassie e degli ammassi di galassie
X-ray image of the
Virgo Cluster
X-ray image of the
Coma Cluster
16 Mpc* faraway
100 Mpc faraway
2.5° radius (compare to
moon 0,25°)
More than 1000 galaxies
*1 Mpc = 1019 km
X-ray image of a
Simulated Cluster
Standard CDM model
2563 simulation
Temperature of
infalling matter
Temperatures range from 0
to 108 K
…e così via…
Firenze, 10-11 Giugno 2003, C. Gheller
Quanto costa una simulazione ???
le strutture cosmologiche sono influenzate da processi che si sviluppano su
un range di scale enorme; alcuni esempi:
•le stelle si formano su scale ~ 1 parsec (1016km)
•le galassie sono su scale ~10 kparsec
•le forze di pressione influenzano scale ~ 10 Mparsec
•la gravità agisce sino a scale comparabili con quella dell'universo, ~ 1000
Mparsec
Il range dinamico ideale sarebbe quindi di 109 !!!!!!!!
Vediamo cosa significa questo in termini di memoria del calcolatore:
1283
Min. 170 Mbyte
Qualche ora/giorno
PC
5123
Min. 10 Gbyte
Qualche settimana
Server multiprocessore
10243
Min. 80 Gbyte
Qualche mese
Supercomputer
Firenze, 10-11 Giugno 2003, C. Gheller
Descrizione del problema
Le osservazioni ci dicono che l'universo e' costituito da 3 componenti fondamentali:
la materia barionica:
protoni, neutroni - costituisce tutto cio' che emette "luce" ;
la radiazione:
prodotta dai corpi celesti e residuo del Big Bang (radiazione di fondo) - la sua
densita' di energia e' molto minore di quella della materia;
la materia oscura:
natura ignota (neutrini massivi, altre particelle elmentari, pianeti ???...) - scoperta
nelle galassie confrontando le curve di rotazione delle stelle (indice della
distribuzione di massa totale) con le curve di luminosita' (indice della distribuzione di
barioni). Confermata poi da innumerevoli altre osservazioni
Firenze, 10-11 Giugno 2003, C. Gheller
Soluzione del Problema
Focalizziamoci sulle due componenti materiali:
la materia oscura ha solo interazioni a lungo raggio (gravita'). Per essa e' adeguato
un approccio numerico tipo N-Body:
la materia e' rappresentata come un insieme di particelle il cui comportamento e'
descritto
dalla
soluzione
delle
relative
equazioni
del
moto
la materia barionica risente anche di interazioni a corto raggio:
essa quindi puo' venire piu' adeguatamente descritta come un fluido, la cui
evoluzione viene caratterizzata dalla soluzione di un sistema di equazioni che
stabiliscono la conservazione di
•materia
•quantita' di moto
•energia
Infine e' necessario calcolare il campo gravitazionale che accoppia le due
componenti.
Firenze, 10-11 Giugno 2003, C. Gheller
Equazioni del modello
Materia oscura: equazioni del moto.
Materia Barionica: equazioni di conservazione.
+ equazione di stato
Campo e Forze Gravitazionali: equazione di Poisson.
Firenze, 10-11 Giugno 2003, C. Gheller
L’algoritmo Fluidodinamico alle Differenze Finite
Con l'approccio alle differenze finite spazio e tempo vengono rappresentati in modo
discreto, come griglia computazionale multidimensionale.
Su questa griglia le vengono definite le quantità caratteristiche del sistema. Queste
possono venire interpretate come:
•valori medi sulle cella della griglia,
•valori puntuali sui nodi della griglia.
Sulla griglia spazio-temporale, vengono riscritte in modo approssimato le equazioni
che caratterizzano la dinamica del sistema, in modo che esse dipendano dai valori
delle quantità caratteristiche approssimati sulla griglia stessa.
Firenze, 10-11 Giugno 2003, C. Gheller
L’algoritmo Fluidodinamico alle Differenze Finite
Con l'approccio alle differenze finite spazio e tempo vengono rappresentati in modo
discreto, come griglia computazionale multidimensionale.
L'idea base del modello numerico e' la seguente:
La materia barionica viene descritta come un fluido;
•il fluido e' descritto attraverso la sua densita', la sua pressione, la sua velocita';
•lo spazio viene rappresentato da una griglia di NxNxN elementi;
•le variabili del fluido vengono calcolate sulle celle della griglia;
•i valori calcolati sono i valori medi della variabili sulle celle;
•questi valori sono calcolati risolvendo con una qualche approssimazione numerica
le equazioni di conservazione:
...ad esempio, al primo ordine
(scarsa accuratezza),
in una dimensione, l'equazione di
conservazione della massa si può
approssimare alle differenze finite
come:
Firenze, 10-11 Giugno 2003, C. Gheller
Parallelizzazione di algoritmi fluidodinamici alle differenze finite
A(i , j)
E’ un problema “semplice”.
Per loro intrinseca natura le equazioni fluidodinamiche
sono locali.
Pertanto la parallelizzazione del codice si può
ottenere in modo efficiente attraverso una
suddivisione a blocchi del dominio, tipo quella
presentata per l’analisi delle immagini.
Il codice risultante scala molto bene ed e’ molto ben
bilanciato.
Algoritmi di questo tipo si prestano molto bene alla
parallelizzazione !!!
Firenze, 10-11 Giugno 2003, C. Gheller
Esempio:
Data
la
località
della
soluzione il processo N ha
bisogno di conoscere SOLO
la prima riga di dati dei
processi N-1 e N+1
Si può dividere il dominio
computazionale a “bande”.
PE N+1
Ghost Region
PE N
Ghost Region
PE N-1
L’array contenente la variabile viene definito con 2 righe in più,
contenenti le “Ghost Regions”. La comunicazione coinvolge
esclusivamente le ghost regions e i bordi. Il resto del calcolo è
completamente locale.
Firenze, 10-11 Giugno 2003, C. Gheller
L’algoritmo N-Body
Nel caso di simulazioni cosmologiche di sola materia oscura, si considera la sola
massa, quindi le sole interazioni gravitazionali. In generale comunque la
procedura puo' essere generalizzata a casi più articolati e complessi:
•ogni particella ha massa M data;
•la velocità di ogni particella viene calcolata in istanti di tempo successivi
(discretizzazione del tempo) risolvendo le equazioni del moto con
opportune tecniche alle differenze finite.
•La forza gravitazionale che agisce su ciascuna particella viene calcolata
attraverso un'opportuna approssimazione numerica.
Firenze, 10-11 Giugno 2003, C. Gheller
Dinamica delle Particelle
Equazioni del moto
Approx. del primo ordine
Approx. del secondo ordine: Leapfrog
Approx. del secondo ordine:
Lax Wendroff (two steps)
Firenze, 10-11 Giugno 2003, C. Gheller
Calcolo delle Forze
Approx. del primo ordine
Metodo Particle Particle (PP):
la forza su una particella è la risultante di tutte le interazioni con le altre particelle
Vantaggi - accuratezza, risoluzione
Limiti - numero di particelle (scala con N2)
Metodi Tree:
la forza è lo sviluppo di multipolo dovuto alla
distribuzione delle particelle
Vantaggi - risoluzione, scalabilità (N log(N))
Limiti - complessità della costruzione dell'albero
Firenze, 10-11 Giugno 2003, C. Gheller
Metodo Particle Mesh
la forza viene calcolata attraverso una doppia trasformata di Fourier
sul campo di densità associato alla distribuzione di particelle ...
Vantaggi - scalabilità, rapidità, semplicità dell'algoritmo
Limiti - Risoluzione spaziale
1)
2)
3)
4)
discretizzo lo spazio con una griglia cubica si
lato L, numero di punti per dimensione N e
(quindi) risoluzione spaziale dx=L/N.
calcolo il campo di densità associato alla
distribuzione di particelle distribuendo la massa
di ciascuna particella sui punti griglia.
risolvo l'equazione di Poisson (che lega la
densità al potenziale gravitazionale) attraverso
una procedura FFT. In questo modo calcolo il
campo gravitazionale sulla griglia.
calcolo la forza differenziando il potenziale e
stimo la forza che agisce sulla particella con
una procedura analoga alla 2).
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione dei dati
Distibuzione dei dati tra i processori:
In questo caso esistono due tipi di vettori
•quelli che descrivono le particelle (posizione, velocita');
•quelli che descrivono i campi sulla griglia (potenziale gravitazionale, densita').
questi ultimi vengono distribuiti a piani paralleli. Le particelle vengono suddivise tra i
processori:
NB: non c'e' correlazione tra la posizione delle
particelle locali ad un processore, l'id del
processore, e il piano della griglia locale a quel
processore
Questo è un problema essenziale per la parallelizzazione !!!
Firenze, 10-11 Giugno 2003, C. Gheller
Firenze, 10-11 Giugno 2003, C. Gheller
Firenze, 10-11 Giugno 2003, C. Gheller
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh:
considerazioni.
1)
in generale una particella sarà locale ad un processore diverso da quello che
possiede i dati su griglia necessari a calcolarne la dinamica.
2)
il numero di particelle che "migrano" cambia continuamente e in modo
disomogeneo man mano che la simulazione procede.
questo porta ad un algoritmo che in generale richiede molta comunicazione e forti
sbilanciamenti nel carico di lavoro che ogni processore deve realizzare !!!
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro.
La scelta delle modalità di distribuzione del lavoro è quindi particolarmente delicata
per l'efficienza del codice.
Nel caso più semplice, si può decidere se far svolgere il lavoro al processore che ha
locale la particella o quello che ha locali i dati su griglia.
Se lo svolge il PE che "possiede" la particella devo fargli avere tutti i dati delle celle
vicine: gli algoritmi standard di distribuzione della massa o di calcolo delle forze
richiedono da 8 a 27 celle.
In questo modo comunque ogni processore calcola la dinamica di NPART_TOT/NPES
particelle (perfetto bilanciamento).
Se il lavoro lo svolge il PE che "possiede" le celle questo deve richiedere soltanto 6
informazioni (la posizione e la velocita') della particella.
Tuttavia si potrebbe verificare che un solo PE esegue tutto il lavoro mentre gli altri
attendono inoperosi.
La scelta quindi dipende dal bilancio tra queste due opposte esigenze: assicurare
che la comunicazione sia minima e garantire un buon bilanciamento del carico di
lavoro.
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Nel esempio che segue, abbiamo invece implementato la seconda strada, che
può portare a forti sbilanciamenti, ma ottimizza la comunicazione.
L'implementazione è basata su MPI. Consta dei seguenti step:
•Individuazione della posizione della particella sulla griglia
•Creazione di liste di particelle associate ai vari processori
•Comunicazione dei dati delle particelle tra i vari processori
•Calcolo della densita'
•Calcolo del potenziale
•Calcolo della forza
•Integrazione delle equazioni del moto
La strategia è quindi quella di sapere in anticipo quali particelle cadono in
processori diversi da quelli a cui appartengono e comunicare le relative
informazioni, in un unico processo di comunicazione, ai processori di
competenza che a questo punto possono eseguire il calcolo in modo
completamente locale
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Questo tipo di implementazione consente il controllo completo della
comunicazione ma necessita di sostanziali modifiche rispetto alla versione
sequenziale.
Step 1: creo alcuni array per lo scambio di informazioni “generali” tra i
processori
allocate (count_pe_s(npes),STAT=error)
allocate (count_pe_r(npes),STAT=error)
count_pe_s=0
count_pe_r=0
Ogni posizione in questi array
contiene il numero di particelle
che
dovranno
essere
inviate/ricevute
a/da
un
processore remoto.
!
! find how many particles fall in each processor
!
n_tot_s=0
do ipart=1,nparmax
In questo loop l’array viene
i3=int(x3(ipart)+0.5)
effettivamente riempito.
if(i3.eq.0)i3=npes*nz
ipe3=(i3-1)/nz
if(ipe3.ne.mype)then
count_pe_s(ipe3+1)=count_pe_s(ipe3+1)+1
n_tot_s=n_tot_s+1 Conto quante particelle in
endif
totale ogni processore invia
enddo
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Step 2: comunico il numero di particelle che i vari processori si scambiano.
CALL MPI_Alltoall ( count_pe_s, 1, MPI_INTEGER, count_pe_r, 1,
MPI_INTEGER, MPI_COMM_WORLD, ierr )
…or, using send end recv…
do i=1,npes
if(i-1.ne.mype)then
to_pe=i-1
from_pe=i-1
CALL MPI_ISend(count_pe_s(i),1,,to_pe,10,&
MPI_COMM_WORLD,ierr)
CALL MPI_IRecv( (i),1,,&
from_pe,10,,status,ierr)
endif
enddo
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Step 3: Sapendo quante particelle devono venire comunicate, preparo i buffer
con cui comunicare le rispettive posizioni e velocità.
allocate
allocate
allocate
allocate
allocate
allocate
(x1_aux_s(n_tot_s+1),STAT=error)
(x2_aux_s(n_tot_s+1),STAT=error)
(x3_aux_s(n_tot_s+1),STAT=error)
(v1_aux_s(n_tot_s+1),STAT=error)
(v2_aux_s(n_tot_s+1),STAT=error)
(v3_aux_s(n_tot_s+1),STAT=error)
do ipart=1,nparmax
i3=int(x3(ipart)+0.5)
if(i3.eq.0)i3=npes*nz
ipe3=(i3-1)/nz
if(ipe3.ne.mype)then
index0=point_pe_aux(ipe3+1)
x1_aux_s(index0)=x1(ipart)
x2_aux_s(index0)=x2(ipart)
Step 4: Riempio i buffer
x3_aux_s(index0)=x3(ipart)
v1_aux_s(index0)=v1(ipart)
v2_aux_s(index0)=v2(ipart)
v3_aux_s(index0)=v3(ipart)
n_pepe(index0)=ipart
point_pe_aux(ipe3+1)=point_pe_aux(ipe3+1)+1
endif
enddo
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Step 5: Allo stesso modo alloco i vettori che dovranno ricevere i dati.
allocate
allocate
allocate
allocate
allocate
allocate
(x1_aux_r(n_tot+1),STAT=error)
(x2_aux_r(n_tot+1),STAT=error)
(x3_aux_r(n_tot+1),STAT=error)
(v1_aux_r(n_tot+1),STAT=error)
(v2_aux_r(n_tot+1),STAT=error)
(v3_aux_r(n_tot+1),STAT=error)
do i=1,npes
if((i-1).ne.mype)then
to_pe=i-1
from_pe=i-1
index0=point_pe(i)
index1=point_pe_r(i)
stride0=count_pe_s(i)
stride1=count_pe_r(i)
…
CALL MPI_Irecv(x1_aux_r(index1),stride1,MPI_DOUBLE_PRECISION,&
from_pe,20,MPI_COMM_WORLD,req(1),ierr)
CALL MPI_Isend(x1_aux_s(index0),stride0,MPI_DOUBLE_PRECISION,&
to_pe,20,MPI_COMM_WORLD,req(2),ierr)
…
CALL MPI_WAITALL(12,req,status_array,ierr)
endif
enddo
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Step 6: gli array di send non servono più. Per risparmiare spazio li elimino…
deallocate
deallocate
deallocate
deallocate
deallocate
deallocate
(x1_aux_s)
(x2_aux_s)
(x3_aux_s)
(v1_aux_s)
(v2_aux_s)
(v3_aux_s)
Step 7: tutti i dati sono ora locali ai processori. Posso procedere con il calcolo
come se fosse sequenziale.
Calcolo la densità a partire dalle
call density(x1_aux_r,x2_aux_r,x3_aux_r,n_tot) particelle
call fields
Calcolo il campo gravitazionale
call gravity
Calcolo le forze gravitazionali
call nbody(n_tot,x1_aux_r,x2_aux_r,x3_aux_r,&
v1_aux_r,v2_aux_r,v3_aux_r)
Sposto le particelle
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI
Step 8: ripeto al contrario quanto fatto in precedenza, riportando le informazioni
su posizioni e velocità delle particelle comunicate ai processori originari.
Step 9: dealloco tutti gli array ausiliari necessari alla comunicazione. Questa
operazione è utile a risparmiare memoria ma, soprattutto, è necessaria in
quanto il numero di particelle da comunicare e quindi la dimensione degli
array, cambiano ad ogni step evolutivo.
Firenze, 10-11 Giugno 2003, C. Gheller
Particle Mesh: distribuzione del lavoro. Implementazione MPI 2
MPI 1 non consentendo accessi asincroni (one-side communication) alla
memoria, necessita, in questo caso, di una “infrastruttura” piuttosto complessa.
Con MPI 2 si può invece ottenere un’implementazione molto più snella, simile a
quella Open MP - shared memory vista in precedenza.
call MPI_FENCE(…)
call MPI_WIN_CREATE(phi,8*N**3,MPI_DOUBLE_PRECISION,…)
call MPI_FENCE(…)
do j=1,npart
ip=nint(x1(j))
trovo la posizione della particella sulla griglia
jp=nint(x2(j))
kp=nint(x3(j))
If(la particella è “locale”)then
forza_media=...
else
call MPI_LOCK(MPI_LOCK_SHARED,…)
call MPI_GET(dati necessari per il calcolo della forza)
call MPI_UNLOCK(…)
forza_media=…
endif
v_new = v_old + forza_media * dt aggiorno la velocita'
x_new = x+old + v_new *dt
aggiorno la posizione
enddo
call MPI_WIN_FREE(…)
Firenze, 10-11 Giugno 2003, C. Gheller
Tree Code
In questo caso NON viene introdotta alcuna griglia computazionale...
La forza sulla particella i-esima viene approssimata come somma dei contributi
delle particelle "vicine":
... e di quell lontane, attraverso un'approssimazione di multipolo del potenziale
gravitazionale:
ove con il pedice cm viene indicato il centro di massa dell'intero set di
particelle e con M la massa totale.
Tipicamente si blocca lo sviluppo in serie all'ordine di quadrupolo
Firenze, 10-11 Giugno 2003, C. Gheller
Tree Code
Quali i vantaggi ???
Il multipolo viene calcolato non sulle particelle, ma su un sottoinsieme di celle.
Le celle vengono calcolate secondo una procedura ad albero:
step 1: definisco un angolo critico A
step 2: divido il box computazionale in 8 celle ciascuna di dimensione l. Calcolo
il centro di massa di ciascuna cella e calcolo la distanza d tra la particella e il
centro di massa.
Se l/d > A la cella è lontana e di essa si calcola il quadrupolo,
se l/d < A suddivido ulteriormente la cella
step 3: ripeto la procedura fino a quando rimango con cellette vuote o
contenenti un'unica particella. A queste applico la legge di Newton.
Firenze, 10-11 Giugno 2003, C. Gheller
Tree Code: implementazione parallela
Due tipi di vettori fondamentali:
•Particelle (posizione, velocità, e massa)
•Tree (posizione, massa del baricentro)
La distribuzione dei dati è fatta in modo
da garantire una distribuzione
omogenea della memoria.
Il lavoro è distribuito equamente tra le
particelle (load balancing)
Viene attribuito lo stesso numero di
particelle ad ogni processore.
Utilizzando l'albero le particelle
vengono redistribuite dinamicamente
tra i processori durante l'evoluzione in
modo che particelle fisicamente
vicine siano anche nello stesso
processore)
Le celle dell'albero vengono numerate
progressivamente dalla più grande
alle più piccole.
I livelli più grossolani vengono distribuiti
in modo ciclico, in modo da evitare il
problema di accesso a risorse critiche:
se i livelli più grossolani fossero tutti sullo
stesso processore tutti andrebbero ad
accedere alla memoria di quel
processore.
I livelli più fini sono invece locali al
processore che possiede il "padre".
Firenze, 10-11 Giugno 2003, C. Gheller
Tree Code: implementazione parallela
Problema fondamentale del tree code.
Calcolo della forza: ogni particella deve recuperare i dati remoti di tutte le
particelle "vicine" che non stanno nello stesso processore tutti i quadrupoli
delle celle "lontane"
Calcolo del tree: per calcolare la componenente di quadrupolo della
cella può essere necessario accedere a particelle "remote" (nel senso della
memoria)
La comunicazione è estremamente pesante
vengono utilizzate one side comunication (latenza minima) basate su:
shmem T3E
Origin 3000
LAPI
SP3
Implementazione MPI 2 è in corso
Firenze, 10-11 Giugno 2003, C. Gheller
Tree Code: ottimizzazione
1 Grouping
Particelle molto vicine possono risentire approssimativamente della stessa
interazione di quadrupolo da celle lontane.
Pertanto è inutile ricalcolarla per ogniuna di queste particelle.
E' possibile quindi:
•definire due distanze caratteristiche R e D;
•partendo da posizioni opportune (zone particolarmente dense)
raggruppare le particelle che stanno in una sfera di raggio R;
•calcolare la forza su queste particelle come:
Ftot = Fnear + Ffar, ove:
Fnear viene calcolata in modo "standard" per particelle e celle che distano
meno di D,
Ffar viene calcolata una volta per tutte per particelle che distano più di D.
Firenze, 10-11 Giugno 2003, C. Gheller
Tree Code: ottimizzazione
2 Dynamic Load Balancing
A causa dei diversi tempi di calcolo delle forze per ciascuna particella, un
processore può terminare prima degli altri.
In questo caso (grazie alle one-side comunications, che non presumono
comunicazione) può "aiutare" un'altro processore a completare il suo
carico di lavoro.
3 Data Buffering
Spesso gli stessi dati, relativi a particelle o celle remote, servono per il
calcolo delle forza che agisce su particelle diverse. Pertanto è stato
implementato un meccanismo di "cache". La memoria del processore non
allocata, viene riempita (per quanto possibile) con dati prelevati in
precedenza su altri processori, che a quel punto divengono locali.
Firenze, 10-11 Giugno 2003, C. Gheller