Views

Biuletyn nr 27

From KDMwiki

Jump to: navigation, search
Biuletyn KDM
1 | 2 | 3 | 4 | 5
6 | 7 | 8 | 9 | 10
11 | 12 | 13 | 14
15 | 16 | 17 | 18
19 | 20 | 21 | 22
23 | 24 | 25 | 26
27 | 28 | 29 | 30
31 | 32
Lista biuletynów

Biuletyn nr 27 (12 marca 2008)

Spis treści

Wydarzenia: Warsztaty KFnrD

Autor: Łukasz Bolikowski

W dniach 11-15 lutego 2008 odbyły się w ICM warsztaty modelowania matematycznego i komputerowego dla stypendystów Krajowego Funduszu na rzecz Dzieci. Jednym z celów tego stowarzyszenia jest: "poprawa warunków rozwoju dzieci i młodzieży, zwiększenie szans pełnego rozwoju dzieciom wybitnie uzdolnionym".

Od ponad 10 lat ICM współpracuje z KFnrD, m.in. organizując warsztaty dla stypendystów programu (relacje z ostatnich lat: 2006, 2007). Oprócz ogromnej przyjemności płynącej z pracy z uzdolnioną, ciekawą świata młodzieżą jest także korzyść bardzo praktyczna: młodzi ludzie często uczestniczą w pracach nad rzeczywistymi problemami badawczymi i w niektórych przypadkach - już jako studenci - kontynuują współpracę z ICM. Sześcioro byłych stypendystów Funduszu jest obecnie pracownikami naszego centrum.

W czasie tegorocznych warsztatów uczestnicy mieli do wyboru jeden z pięciu problemów:

  • Diffusion-Limited Aggregation (autor: Michał Łopuszyński)
  • Obliczenia na kartach graficznych (autorzy: Łukasz Ligowski i Jacek Migdał)
  • Analiza obrazów ultrasonograficznych (autor: dr Krzysztof Nowiński)
  • Analiza struktury linków w Wikipedii (autor: Łukasz Bolikowski)
  • Budowa modelu komunikacji krajowej (autor: Franciszek Rakowski).

Uczestnicy wysłuchali ponadto trzech wykładów: prof. Marek Niezgódka opowiedział o modelowaniu matematycznym, dr Anna Trykozko przedstawiła Metodę Elementu Skończonego, a dr hab. Wojciech Grochala pokazał na czym polega praca chemika. Tradycją wszystkich funduszowych warsztatów dotyczących nauk poznawczych są wydarzenia kulturalne: w tym roku Fundusz zaprosił uczestników i organizatorów warsztatów w ICM na "Jezioro Łabędzie" Piotra Czajkowskiego do Teatru Narodowego.

Poniżej znajdują się wybrane wyniki pracy:

Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/9/9d/KFnrD2008_DLA.jpg': No such file or directory.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/9/9d/KFnrD2008_DLA.jpg/400px-KFnrD2008_DLA.jpg'.
Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/e/e0/KFnrD2008_Mapa_przed.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/e/e0/KFnrD2008_Mapa_przed.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/e/e0/KFnrD2008_Mapa_przed.png/400px-KFnrD2008_Mapa_przed.png'.


Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/6/61/KFnrD2008_Mapa_po.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/6/61/KFnrD2008_Mapa_po.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/6/61/KFnrD2008_Mapa_po.png/400px-KFnrD2008_Mapa_po.png'.
Agregacja limitowana dyfuzją (autor: Paweł Lipski) Rozpoznawanie rodzaju terenu przy użyciu sieci neuronowych (autor: Michał Kurtys)
Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/5/50/KFnrD2008_Gestosc_zaludnienia.jpg': No such file or directory.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/5/50/KFnrD2008_Gestosc_zaludnienia.jpg/300px-KFnrD2008_Gestosc_zaludnienia.jpg'.
Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/3/3e/KFnrD2008_Miasta_i_obszary_oddzialywania.jpg': No such file or directory.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/3/3e/KFnrD2008_Miasta_i_obszary_oddzialywania.jpg/300px-KFnrD2008_Miasta_i_obszary_oddzialywania.jpg'.
Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/4/44/KFnrD2008_Natezenie_ruchu.jpg': No such file or directory.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/4/44/KFnrD2008_Natezenie_ruchu.jpg/300px-KFnrD2008_Natezenie_ruchu.jpg'.
Rozpoznawanie miast, obszarów oddziaływania i dróg, oraz model komunikacji krajowej, na podstawie mapy gęstości zaludnienia (autorzy: Jarosław Błasiok i Łukasz Wołochowski)

Ponadto Wojciech Baranowski i Robert Obryk odnaleźli i poprawili niektóre błędne odnośniki w Wikipedii na podstawie analizy podobieństwa struktury artykułów, a Eryk Wieliczko i Paweł Żuk zbliżyli się do teoretycznej wydajności karty graficznej GeForce 8800 podczas analizy podobieństwa tekstów publikacji naukowych.

Superuser: Równoległy rachunek magnetosfer

Autor: Ewa Santon


Jest pożyteczną przyjemnością obserwować fragment świata w myślach lub/i na ekranie. Przyjrzyjmy się tej drugiej możliwości. Fragmentem świata będzie magnetosfera planety, reagująca na płynącą ku niej od Słońca plazmę – zbiór jonów i elektronów. Obraz na ekranie jest wizualizacją wielkości fizycznych, które ten układ opisują. Ich lista jest długa. Na przykład mogą to być: tory wybranych elektronów, ich energia, rozkład w przestrzeni ich gęstości liczbowej, linie indukcji magnetycznej, czy rozkład w przestrzeni natężenia pól – elektrycznego i magnetycznego.


Zmiana stanu magnetosfery - prędkości (v), położeń (r), gęstości prądu (j) i pól – elektrycznego (E) i magnetycznego (B), odbywa się zgodnie z zasadami zachowania, które zapisać można w postaci układu równań różniczkowych:

równania ruchu:

Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \frac{dv}{dt} = \frac{Q}{M}(E + v \times B),


Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \frac{dr}{dt} = v


i równań Maxwell'a:

Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \frac{\delta E}{\delta t} = c^{2} \nabla \times B - \frac{j}{\epsilon_{0}}


Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \frac{\delta B}{\delta t} = - \nabla \times E


Zauważmy, że wartość ilorazu ładunku cząstki (Q) i jej masy (M), Q/M, który występuje w równaniu ruchu jest taka sama dla pojedynczego elektronu, jak dla N-elementowej populacji elektronów, zawartych w małej objętości: Q/M = (N*Q)/(N*M). Zatem takie małe populacje będą poruszać się tak samo, jak pojedynczy elektron. Ten argument uzasadnia zastąpienie wielkiej liczby cząstek nieco mniejszą liczbą tych małych populacji nazwanych makrocząstkami.


Natomiast równania Maxwell'a rozwiązuje się numerycznie na siatce węzłów, a wkład makrocząstek do gęstości prądu w węźle (i, j, k) jest proporcjonalny do części wspólnej makrocząstki i otoczenia węzła, nazywanego komórką. To postępowanie nazywa się metodą cząstki w komórce (ang. Particle-In-Cell).

I tak, po zadaniu początkowych położeń i prędkości cząstek i początkowych rozkładów pól, układ równań jest rozwiązywany w pętli czasowej: w każdym m-tym kroku obliczane są prędkości i położenia cząstek, gęstość prądu i składowe pól. Pomijam problem warunków brzegowych, ale nie mogę nie wspomnieć o ograniczeniach nakładanych przez warunki stabilności numerycznej, szczególnie ten, dotyczący równania ruchu. Z warunku tego wynika, że krok czasowy powinien być na tyle mały, aby cyklotronowa składowa ruchu elektronów była koniecznie obecna w rachunku.

Oto przykład takich obliczeń. Jeśli decydujemy się na to, aby obserwować otoczenie planety o wymiarach 25.6x12.8x12.8 (jednostką długości jest promień planety Rp), z rozdzielczością przestrzenną 0.1Rp, to w układzie jest prawie 4.2 mln węzłów. Jeśli w każdej komórce, poza wnętrzem planety, umieścimy parę {elektron, proton}, to w układzie jest prawie 8 mln cząstek. Jeśli wartość maksymalna indukcji magnetycznej w układzie Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): B_{0} = 1 nT , to krok czasowy nie może być większy niż Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): 2./(Q/M)B_{0} = 0.0112s . Przyjmijmy Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \Delta t = 0.01s . Mamy zatem sześć wektorów dla położeń i prędkości, każdy o rozmiarze 8 mln, oraz sześć tablic dla składowych pól, każda o rozmiarze 4.2 mln. Aby prześledzić 3-minutową ewolucję układu, trzeba wykonać 18 tys. kroków. Obliczenia powinny być wykonane w rozsądnym czasie. Dlatego program został zrównoleglony i był liczony na 'tornado'.

Oto kilka obrazów ilustrujących wybrane momenty ewolucji magnetosfery. Zacznijmy od obszaru symulowanego.


Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/7/7f/Santon01.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/7/7f/Santon01.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/7/7f/Santon01.png/400px-Santon01.png'.
Rys. 1. Domena symulacji; warstwa elektronów, o grubości 2Rp , t = 102s.


Od strony Słońca przez ściankę x = 0 płyną do wnętrza domeny elektrony i protony wiatru słonecznego. Środek planety znajduje się w punkcie [6.3,6.3,6.3]Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): R_{p} . Na rysunku zaznaczono półkulę planety od strony zmierzchu, linie Słońce-planeta i południe-północ oraz ślad dwóch sond. Jedna krąży wokół planety na wysokości Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): 0.25R_{p}

, druga przecięła stronę nocną magnetosfery. Energie elektronów zakodowano kolorem. Maksymalna możliwa wartość energii to 800 eV i tej energii odpowiada  kolor czerwony. Początkowa energia jest na poziomie 1 eV.   

Wcześniej, w chwili t = 54s, na wieczornych obrzeżach magnetosfery, rozpoczął się proces falowy – oscylacje składowych pola magnetycznego. Trzy kolejne obrazy ilustrują jego przebieg.


Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/0/01/Santon02.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/0/01/Santon02.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/0/01/Santon02.png/400px-Santon02.png'.
Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/5/51/Santon03.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/5/51/Santon03.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/5/51/Santon03.png/400px-Santon03.png'.
Rys. 2. Proces falowy na obrzeżach magnetosfery; chwile t = 51s i 81s. Natężenie składowej z-owej pola B w płaszczyźnie Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): z = 6.3R){p} . Środek i kontur planety zaznaczono kropkami.


Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/1/15/Santon04.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/1/15/Santon04.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/1/15/Santon04.png/400px-Santon04.png'.
Rys. 3. Zapis z-owej składowej indukcji magnetycznej B w punkcie [8.1, 4.3,6.3]Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): R_{p} (linia niebieska) oraz, odpowiednio przeskalowane, 1-sza - 0.13 Hz (czerwona) i druga - 0.26 Hz (zielona), składowe gęstości spektralnej tego sygnału.


Łatwo dostrzec oscylacje o okresie 3.5 – 8.5s, a to oznacza obecność częstości f odpowiednio w zakresie 0.28 – 0.11 Hz. Odnotujmy wszakże: w obszarach, gdzie B = 0.1 – 0.2 nT, częstość cyklotronowa elektronów, fce = 2.8 – 5.6 Hz. Zatem f ≈ 0.1 Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): f_{ce}

. Tego typu aktywność falowa jest typowa dla obrzeży magnetosfer.

Jeszcze jeden szczegół z dynamiki magnetosfer. Na Rys. 1 naniesiono także populację elektronów z chwili t = 128s. Energia tych cząstek jest większa od 600 eV. Weszły one do układu w chwili t = 81s, wraz z międzyplanetarnym polem elektrycznym. Są nieliczne, ale rejestrowane przez wirtualną sondę; zapis gęstości liczbowej tych wtórnych elektronów umieszczono na Rys. 4.


Error creating thumbnail: convert: unable to open image `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/2/20/Santon05.png': No such file or directory.
convert: unable to open file `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/2/20/Santon05.png'.
convert: missing an image filename `/media/wwwold.icm.edu.pl/kdmwiki-dir2/images/thumb/2/20/Santon05.png/400px-Santon05.png'.
Rys. 4. Gęstość liczbowa elektronów zarejestrowana na pokładzie sondy wirtualnej. Długość trajektorii znormalizowano do jedności. Symbolem CA oznaczono punkt o najmniejszej odległości od powierzchni planety (ang. closest approach)


Takie zapisy są znane z pomiarów zarejestrowanych przez analizatory cząstek, umieszczane na pokładach sond próbkujących magnetosfery planet. Podsumowując te, z konieczności krótkie, rozważania, można powiedzieć, że istnieją procesy magnetosferyczne, które odtwarza sformułowany na początku model, a to stanowi o jego przydatności. Aby z niego korzystać, potrzebny jest jedynie wielordzeniowy komputer dużej mocy.

Octave - czyli darmowy Matlab

Autor: Maciej Szpindler

Octave to darmowy klon znanego pakietu do obliczeń matematycznych MATLAB. Program jest udostępniany na zasadach licencji GNU i można znaleźć go w większości popularnych dystrybucji linuxa np. Fedora, SUSE, Ubuntu. Filozofia i funkcjonalność programu Octave jest bardzo podobna do tej znanej z MATLABa (zob. artykuł o MATLAB w Biuletynie nr 7 ). Jest to więc środowisko wraz z własnym językiem programowania wysokiego poziomu zaprojektowanym do łatwego wykonywania obliczeń, przede wszystkim algebry liniowej, ale nie tylko. Jest to bardzo wygodny język - pozwala na szybkie zapisanie i przetestowanie własnego programu.

Octave został stworzony tak, aby zapewnić możliwie dużą zgodność ze swoim pierwowzorem - MATLABem. I tak rzeczywiście jest; oczywiscie niektórych zaawansowanych funkcji dostępnych w pakietach rozszerzających - tzw. toolboxach - MATLABa brakuje w Octave, ale można ich szukać w Octave-forge. To zbiór wielu dodatków i usprawnień do Octave rozwijanych trochę "z boku" od głównego nurtu rozwoju programu.

Zatem programy i funkcje zapisane w języku Octave umieszczamy w plikach .m (jak MATLAB). Podamy tu kilka krótkich przykładów jak napisać własny program z użyciem Octave.

Przykłady użycia Octave

Jako przykład obliczymy przybliżoną wartość liczby Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \pi , korzystając ze wzoru:

Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \pi = 4 \int_0^1 \frac{1}{1+x^2}


i całkowania numerycznego metodą prostokątów. Podobnie jak robiliśmy to już w biuletynie przy okazji omawiania Co-Array Fortran i Unified Parallel C. Przykładowy programik może wyglądać następująco:

function [pi_comp,err] = Pi_seq(N)
  
  width = 1/N; 
  Sum = 0;           
  i = 0 : N-1;
  x = (i + 0.5)*width;
  Sum = sum( 4./(1 + x.^2) );
  Sum = Sum/N ;
  err  = Sum - pi;
  pi_comp   = Sum;

end


Kod programu można znaleźć na komputerze halo w pliku /opt/octave/mpitb/Pi/Pi_seq.m.

Octave równoległy

Do obliczeń równoległych z poziomu Octave służy dodatkowy pakiet MPITB. Pozwala on wywoływać funkcje biblioteki MPI z Octave.

Zaczniemy od najprostszego programu - równoległego Hello world. Program może wyglądać następująco Hello.m:

  info       =  MPI_Init;
 [info rank] =  MPI_Comm_rank (MPI_COMM_WORLD);
 [info size] =  MPI_Comm_size (MPI_COMM_WORLD);
 [info name] =  MPI_Get_processor_name;

 fprintf("Hello, MPI_COMM_world! I'm rank %d/%d (%s)\n", rank, size, name);

 info       =  MPI_Finalize;

 MPITB_Clean;

Możemy przetestować program na klastrze halo, mamy tam dostępny pakiet Octave; Program Hello znajdziemy w katalogu: /opt/octave/mpitb/Hello.

Aby uruchomić program na wielu procesorach musimy umieścić zadanie w systemie kolejkowym, przygotowując odpowiedni skrypt octave-mpitb.pbs:

#!/bin/csh
#PBS -S /bin/csh
#PBS -N octave-mpitb 
#PBS -l cput=00:10:00
#PBS -l mem=10mb
#PBS -l nodes=2:ppn=2
setenv MPITB_HOME /opt/octave/mpitb/
use_openmpi_gnu64 
use_octave
cd /opt/octave/mpitb/Hello
mpiexec -np 4 --mca btl self,sm,tcp octave -q --eval Hello

Skrypt zakłada, że użyjemy 4 procesorów w konfiguracji 2 x 2 procesory w węzłach klastra (więcej: tu) Umieszczamy go w kolejce poleceniem:

qsub octave-mpitb.pbs

Program po wykonaniu powinien zwrócić podobny wynik (w odpowiednim pliku zwróconym przez system kolejkowy):

Hello, MPI_COMM_world! I'm rank 1/4 (n96)
Hello, MPI_COMM_world! I'm rank 3/4 (n97)
Hello, MPI_COMM_world! I'm rank 0/4 (n96)
Hello, MPI_COMM_world! I'm rank 2/4 (n97)

Możemy teraz spróbować zrównoleglić nasz program przybliżający liczbę Parser nie umiał rozpoznać (Missing texvc executable; please see math/README to configure.): \pi

w postaci Pi.m:
function [pi_comp, err] = Pi(N)
  MPI_Init;                  
  [info rank] = MPI_Comm_rank (MPI_COMM_WORLD); 
  [info size] = MPI_Comm_size (MPI_COMM_WORLD);

   SLAVE = logical(rank);                 
   MASTER = ~ SLAVE;                        

 width = 1/N; 
 lsum = 0;            
 i = rank : size : N-1;                
 x = (i + 0.5)*width;              
 lsum = sum( 4./ (1+x.^2) );        

 TAG = 7;
 if SLAVE                              
   MPI_Send(lsum, 0,TAG,MPI_COMM_WORLD);
 else                                
   Sum = lsum;                  
   for slv= 1 : size - 1                   
     MPI_Recv(lsum,MPI_ANY_SOURCE,TAG,MPI_COMM_WORLD);
     Sum += lsum;                  
   end                               
 end

 MPI_Finalize;
 MPITB_Clean;          
 if MASTER
   Sum = Sum/N ;                     
   err  = Sum-pi;
   pi_comp  = Sum;
 end
 
end

Program Pi znajdziemy w katalogu: /opt/octave/mpitb/Pi. Odpowiedni skrypt do systemu kolejkowego może wyglądać następująco:

#!/bin/csh
#PBS -S /bin/csh
#PBS -N octave-mpitb 
#PBS -l cput=00:10:00
#PBS -l mem=10mb
#PBS -l nodes=2:ppn=2
setenv MPITB_HOME /opt/octave/mpitb/
use_openmpi_gnu64 
use_octave
cd /opt/octave/mpitb/Pi
mpiexec -np 4 --mca btl self,sm,tcp octave -q --eval "Pi(1e7)"

Zobacz też