ALGORITMA TREE BERHIRARKI UNTUK KOMPUTASI PARTIKEL FISIKA Slamet Santosa*, Aminus Salam*
ABSTRAK
ALGORITMA TREE BERHIRARKI UNTUK KOMPUTASI PARTIKEL FISIKA. Simulasi partikel N-Bodi telah banyak digunakan untuk mempelajari berbagai sistem fisika seperti fisika plasma, fisika akselerator, sistem dinamika molekuler dan astrofisika yang menggunakan interaksi antar partikel sebagai dasar komputasinya. Walaupun perkembangan teknologi komputer memungkinkan untuk mensimulasi banyak partikel yang berinteraksi secara elektrostatik ataupun gravitasi, tantangan untuk mengembangkan algoritma dari komputasi tersebut masih terbuka. Telah dikembangkan metode tree berhirarki yang lebih baik daripada metode titik jala sehingga orde komputasinya mencapai O (N log N). Pada makalah ini dibahas implementasi yang efisien dari algoritma tree berhirarki menggunakan teknik quick-sorting untuk penpartisi domain partikel pada pembentukan tree dan penjejakan tree dengan metode Breath First Search dalam menghitung gaya pada tiap partikel.
ABSTRACT HIERARCHICAL TREE ALGORITHM FOR PARTICLE PHYSICS COMPUTATION. N-Body simulations have been used to study a wide variety of physical system such as plasma physics, accelerator physics, molecular dynamics and astrophysics where a large number of interacting particles are the basic of such simulations. Although computer technologies have made the particle simulations feasible, this type of N-Body problems including those of electrostatic or gravitational forces still present a challenge. The hierarchical tree method offers the possibility of computing the interaction between N particles in a time O(N log N) without using a grid. The efficient implementation of the algorithm using a divide-and-conquer technique based on a quick-sorting technique in the building tree and the tree traversal based on Breath First Search for force calculations are discussed.
PENDAHULUAN Komputasi partikel fisika merupakan basis dari simulasi komputer di bidang fisika plasma, astrofisika, akselerator dan sistem dinamik molekuler [1, 2]. Komputasi yang dimaksud adalah penyelesaian algoritmik dari lintasan banyak partikel pada *
Pusat Penelitian dan Pengembangan Teknologi Maju - BATAN
sistem fisika klasik dan dikenal sebagai masalah N-Body, di mana diberikan suatu domain yang memuat banyak partikel yang saling berinteraksi secara elektrostatik ataupun medan gravitasi. Persamaan gerak dari lintasan banyak partikel sangat sederhana akan tetapi tidak mempunyai penyelesaian analitik yang efisien sehingga hanya efektif apabila dilakukan komputasi dengan komputer berkecepatan tinggi. Pada sistem dengan N buah partikel yang diketahui masing-masing posisinya, massa dan kecepatan awalnya akan bekerja gaya elektrostatik atau gaya gravitasi sehingga tiap-tiap partikel akan bergerak pada lintasannya. Lintasan partikel-partikel dapat dimodel dengan persamaan diferensial (PD) yang dapat diselesaikan dengan mengintegrasikan untuk tiap selang waktu tertentu. Sehingga pada tiap selang waktu selama simulasi dapat dihitung posisi dan parameter-parameter yang lain. Gaya pada tiap partikel f i akibat adanya n buah partikel disekelilingnya dihitung sebagai berikut: fi =
n
∑f
i, j
untuk i ≠ j
j =1
Akibat adanya gaya yang bekerja pada sistem, apabila posisi awal suatu partikel xi ( t 0 ) = X 0 maka pada saat t posisi partikel tersebut berpindah dan dapat dihitung dengan persamaan: xi' ( t ) = f ( t , xi ( t 0 )) Algoritma untuk menyelesaikan komputasi tersebut, dengan N = jumlah partikel, t = 0 dan dt = selang waktu tertentu adalah sebagai berikut: while t < tsimulasi for i = 1 to N f (i ) =
∑
N j =1
j ! = i , f (i , j )
call (integrator) update (posisi) end for t = t + dt end while Dengan demikian jelaslah bahwa komputasi tersebut berorde O(N2). Untuk mempercepat proses komputasi telah dikembangkan metode PIC (Particles In Cell), di mana dibentuk titik-titik jala pada domain partikel yang kemudian pada setiap titik jala dihitung potensial dari partikel yang dilingkupi sehingga dengan harga-harga potensial pada tiap titik jala dapat dihitung gaya-gaya dari tiap partikel. Metode PIC ini
mempunyai orde komputasi O(N + M log M) dengan M = jumlah titik-titik jala dan N = jumlah partikel pada sistem. Sehingga metode ini mempunyai efisiensi yang tinggi untuk N yang sangat besar dan distribusi partikel serba sama. Algoritma Barnes and Hut (BH) [3] menghitung langsung gaya-gaya pada N buah partikel pada sistem 3D menggunakan struktur data tree berhirarki dan mempercepat komputasi sehingga berorde O(N log N). Teknik yang dipakai pada algoritma ini adalah divide and conquer, yaitu membagi domain partikel menjadi delapan apabila terdapat dua buah partikel berada pada domain yang sama. Sehingga apabila total domain partikel adalah sebagai root dari tree maka masing-masing subdomain adalah subtree. Pembagian domain dilakukan secara rekursif hingga pada tiap-tiap leaf dari tree terdapat sebuah partikel dan kemudian gaya pada tiap partikel dihitung dengan cara menjejaki oct-tree yang terbentuk. Pada makalah ini dibahas implementasi sekuensial dari algoritma BH tree berhirarki yang meliputi pembentukan tree dan teknik divide and conquer menggunakan metode yang berbasis pengurutan cepat (quick-sorting), penjejakan octtree yang terbentuk dengan metode Breath First Search (BFS) dan analisis kompleksitas orde komputasi. Untuk memudahkan pembahasan, dibahas untuk sistem 2D dimana tree yang terbentuk adalah quad-tree.
LANDASAN TEORI
Pada sistem dengan N buah partikel, Xi = <xi, yi, zi>, Vi =
dan mi = <m1, m2, …, m3> akan bekerja gaya-gaya pada tiap partikel, F = m * A sehingga tiap partikel akan bergerak dan berpindah posisi dari X1(t = 0) ke X2 (t ≠ 0) dengan dX dV laju perubahan dan laju kecepatan gerak . Hubungan antara X(t) dan V(t) dt dt menyatakan bahwa sistem adalah merupakan PD rangkap, dan apabila diberikan suatu harga awal X(t0) = X0, penyelesaiannya dapat dituliskan sebagai berikut: X ' ( t ) = f ( t , X (t )),
d ( t , X (t )) = A; X (t 0 ) = X 0 dt
Pada masalah medan gravitasi [5] yang dipergunakan untuk uji masalah dalam implementasi algoritma ini, gaya yang bekerja pada satu partikel mengikuti hukum gravitasi, sehingga percepatan yang diderita oleh partikel tersebut adalah:
Ai =
N X j − Xi G * mj * mi j =1, j ≠i r 3i , j
∑
merupakan fungsi percepatan derivatif yang mempunyai dua buah argumen massa dan posisi. Dengan demikian harga ini harus dihitung sebelum proses update posisi untuk setiap step waktu berikutnya. Untuk menghitung fungsi kecepatan derivatif dipergunakan sebuah integrator. Metode Euler adalah integrator sederhana yang menghitung X n +1 = X n + h(t , X n ) dimana h dapat diperoleh menggunakan deret Taylor. Metode lain adalah metode Runge-Kutta, yang merupakan perluasan metode Euler dan dapat menghitung dengan dX lebih teliti. Metode ini menggunakan step penambahan kecil pada fungsi derivatif dt untuk menghitung harga X n+1 . Standar metode Runge-Kutta menggunakan empat buah step penambahan kecil k1, k2, k3 dan k4 sebagai parameter-parameter yang digunakan untuk menghitung harga X n+1 dengan persamaan:
1 X n +1 = X n + ( k 1 + 2 k 2 + 2 k 3 + k 4 ) 6 Dalam hal persamaan gerak partikel pada program yang dibuat, percepatan d (t , X ( t )) A= adalah fungsi Accel(mj, Xj) dengan mj = massa partikel-partikel lain dt dengan posisi Xj dan harga-harga k1, k2, k3 dan k4 dihitung sebagai berikut: k 1 = h( t n , X n ) = Vn k 2 = h (t n +
k h h , X n + 1 ) = Vn + * Accel ( m j , X j ) 2 2 2
k 3 = h (t n + h , X n +
k2 h h ) = Vn + * Accel ( m j , X j + * Vn ) 2 2 2
h h * Vn + * Accel ( m j , X j )) 2 2 Pada distribusi partikel yang tidak serba sama, akan didapati relatif kumpulan partikel, yang dibenarkan pada saat menghitung gaya, dianggap sebagai partikel tunggal apabila jaraknya relatif jauh dari partikel yang ditinjau [4]. Algoritma BH menggunakan cara tersebut dan dengan teknik open-test faktor ketelitian θ untuk menghitung gaya-gaya pada setiap partikel sebagai berikut. Pada sistem 2D, domain k 4 = h(t n + h, X n + k 3 ) = Vn + h * Accel( m j , X j +
D , D = panjang sel dan r = jarak r dari partikel yang ditinjau ke pusat sel lebih kecil dari faktor ketelitian θ , maka partikel-partikel yang ada di dalam sel adalah partikel tunggal berposisi di pusat sel, dengan massa sama dengan total massanya. Proses pembentukan struktur tree berhirarki dari algoritma BH pada sistem 2D adalah sebagai berikut. Dibangkitkan partikel secara random pada domain partikel yang didefinisikan, apabila pada sektor yang sama terdapat dua buah partikel, sektor tersebut dibagi empat subsektor yang masing-masing sebagai “anak” pada struktur tree. Dengan demikian setiap struktur sektor adalah node dari tree yang terbentuk. Proses tersebut diulang secara rekursif sehingga pada setiap sektor terdapat satu buah partikel. Pada struktur sektor dilengkapi empat buah pointer, sehingga terbentuk struktur tree berhirarki, dimana root adalah domain yang memuat seluruh partikel. Jika sebuah node adalah v, maka *p(v) adalah subset dari node v, sehingga * p(v ) = U ki =1 p( vi ) untuk tiap-tiap k = 2-4. Gambar 1 memperlihatkan divide and conquer dari struktur oct-tree dari algoritma BH untuk sistem 3D. Untuk menghitung gaya pada setiap partikel dilakukan penjejakan struktur tree (tree traversal) dengan mengikuti pointer pada tiap sektor. partikel dibagi menjadi banyak sel, apabila rasio
Gambar 1. Hasil Divide and Conquer pada Sistem 3D
METODOLOGI Pada algoritma BH secara implisit menjelaskan bahwa pembentukan tree dimulai dengan didefinisikan domain partikel berupa kubus besar yang cukup memuat partikel dari sistem. Partikel kemudian dibangkitkan secara acak dan dimasukkan satu persatu, jika ada dua buah partikel menempati bidang yang sama, bidang tersebut dibagi delapan dan diaktifkan root pointer menunjuk delapan subbidang tersebut dan seterusnya. Pada implementasi algoritma BH ini, pembentukan tree dikembangkan dan dimulai dengan mendefinisikan domain partikel dan semua partikel dari sistem dimasukkan ke dalamnya sebelum pembagian domain partikelnya. Hal ini ternyata mempercepat proses pembentukan tree karena dapat dipakai teknik pengurutan cepat.
Pembentukan Tree Didefinisikan P adalah sejumlah partikel dari sistem pada Rd dengan d adalah konstanta yang menyatakan dimensi. Untuk d = 2, sistem adalah 2D, sebuah sektor pada sembarang level h = i adalah produk Cartesian s = [ x i , yi ] ∈ R 2 yang memuat sejumlah partikel. S[P] adalah root dengan level h = 0. Domain root didefinisikan dengan l x [ P ] dan l y [ P] yang masing-masing adalah panjang bidang x dan y. Titik sudut kiri bawah dari tiap-tiap sektor didefinisikan sebagai (xl, yl) sehingga jumlah partikel termuat dalam suatu sektor dapat dicari dengan (xl, yl, xl + lx[P], yl + ly[P]). Teknik divide and conquer membagi tiap sektor melalui titik tengah sektor, karena pada program yang dibuat ditentukan lx[P] = ly[P] maka titik-titik tengah tiap sektor adalah c[ s] = ( x l , y l ) + l[ P ] / 2 . Data partikel disimpan pada sebuah array. Data posisi partikel disamakan dengan domain partikel dengan cara memetakan setiap partikel dengan memproyeksikan data posisi melalui variabel-variabel domain partikel. Proses selanjutnya adalah melakukan pengurutan secara parsial pada data partikel dengan cara membandingkan setiap koordinat partikel dengan titik tengah sektor c[s] dari sektor yang ditinjau. Dengan demikian dapat ditentukan segmen-segmen pemisah pada array yang menentukan subsektor1 sampai dengan subsektor4, yang masing-masing dinyatakan sebagai s[p1], …, s[p4]. Dengan didapatkannya subsektor-subsektor tersebut, masing-masing dapat ditentukan kondisi node-nya dengan cara tes kondisi yaitu dengan menghitung jumlah partikel yang termuat pada subsektor tersebut sebagai berikut: • Jika subsektor tidak memuat partikel maka subsektor adalah node = Nill • Jika subsektor memuat tepat satu partikel maka subsektor adalah node = leaf
•
Jika subsektor memuat banyak partikel maka subsektor adalah node yang padanya harus dilakukan pembagian lebih lanjut.
Struktur node yang dipergunakan pada program adalah sebagai berikut: struct NODE { stuct NODE *subsnodes[4]; int num_subnodes; /* other data */ }; Gambar 2 adalah menyatakan quad-tree yang terbentuk.
Gambar 2. Struktur Quad-tree
Penjejakan Tree Ada beberapa metode yang efisien untuk menjejaki struktur tree untuk mencari node-node pada tree, antara lain metode Depth First Search (DFS) dan medode BFS. Untuk masalah komputasi partikel seperti pada algoritma BH, metode BFS lebih efisien daripada metode DFS karena node-node yang dicari dan memenuhi kriteria faktor ketelitian θ hanya akan dijejaki tepat satu kali, sehingga dengan DFS akan menimbulkan langkah balik yang tidak diinginkan. Prinsip kerja metode BFS untuk menghitung gaya pada tiap partikel dengan faktor ketelitian θ adalah sebagai berikut. Dimulai dengan meninjau root tree, jika root mempunyai subnode yang dinyatakan pada pointernya, ditinjau masing-masing
subnode. Pada tiap-tiap subnode yang ditinjau dilakukan tes kriteria θ , jika dipenuhi maka subnode adalah sebagai partikel tunggal dan penjejakan dihentikan. Sebaliknya, jika kriteria tersebut tidak dipenuhi, penjejakan dilanjutkan pada subnode yang lain. Demikian seterusnya secara rekursif dilakukan penjejakan pada tiap-tiap subnode yang kriteria θ -nya belum dipenuhi sampai dengan dicapai leaf dari tree. Pada penjejakan tree dengan metode BFS untuk menghitung gaya-gaya pada tiap partikel diaktifkan tiga buah List yang fragmen programmnya dapat dijelaskan sebagai berikut: Function BFS((*KriteriaBH)(subnode *)); copy root ke input_work_List; while (!Empty (input_work_List)) { for (each daughter node) { if (KriteriaBH (daughter node)) { copy (daughter node (interaction_List)); else copy (daughter node (output_work_List)); } } input_work_List = output_work_List; }; Total gaya pada suatu partikel dapat dihitung dengan mengevaluasi gaya-gaya dari partikel-partikel pada interaktion_List.
HASIL DAN PEMBAHASAN Proses pembentukan tree pada implementasi algoritma BH ini mempunyai orde komputasi O(N log N) dengan N = jumlah partikel pada sistem, yang dicapai pada kondisi paling buruk (worst-case), yaitu dicapai pada saat distribusi partikel serba sama. Hal tersebut pada sistem 2D, jumlah node yang terbentuk adalah 4 h , dimana h = level dari tree yang terbentuk. Proses penjejakan tree untuk menghitung gaya-gaya pada tiap partikel, untuk harga-harga faktor ketelitian θ yang tidak nol, dapat dianalisis menggunakan
Gambar 3. Distribusi Speris Partikel Berhirarki pendekatan distribusi partikel speris seperti pada gambar 3. Kerapatan distribusi partikel pada sistem dapat dihitung dengan persamaan: 4π 3 n= N /( R ) 3 dimana R adalah radius dari domain partikel yang didefinisikan. Dengan mengorganisasi distribusi partikel di dalam sel-sel berhirarki, maka jumlah interaksi pada satu partikel adalah: nint = n0 +
j .sel
∑ i =0
i n subsel = n0 +
4πri3θ 24 = n0 + 2 4π 3 3 θ ri θ / 8 3
dimana n0 = jumlah interaksi langsung partikel-partikel yang dekat dengan partikel yang ditinjau. Untuk N jumlah partikel pada sistem yang besar, maka jumlah interaksi pada satu partikel dapat dihitung sebagai berikut: nint =
24 log(θ ( 3 N / 4π ) 3 ) 4π 1 N * + ≈ log 2 2 3 log( 1 + θ ) 3 θ θ θ
Dengan demikian, untuk kondisi paling buruk (worst-case) waktu yang diperlukan oleh CPU untuk menghitung gaya pada satu partikel dalam sistem adalah berorde O(log N), yang berarti jumlah operasi untuk menghitung gaya-gaya pada sistem dengan N buah partikel adalah O(N log N). Hasil eksekusi program untuk N = 104 jumlah partikel, dengan berbagai harga faktor ketelitian θ dan dibandingkan dengan algoritma interaksi langsung (PP) diberikan pada gambar 4.
Gambar 4. Eksekusi Program dengan Berbagai Harga θ Dari gambar tersebut dapat dilihat transisi sekala dari orde komputasi O(N2) ke orde komputasi O(N) pada saat harga θ bertambah, yang menyatakan bahwa metode tree berhirarki akan lebih efisien apabila diaplikasikan untuk θ yang diperbesar.
KESIMPULAN DAN SARAN Dengan menggunakan teknik pengurutan cepat pada pembentukan tree dan tiaptiap segmen pembagi dapat digunakan sebagai sentinel pada program, maka eksekusi program lebih efisien dan menjadi lebih cepat. Akan tetapi karena masih menggunakan empat buah pointer pada struktur data node yang dipakai, maka diperlukan jumlah memori yang lebih besar untuk mengalokasi data partikel yang termuat pada tiap node. Hal tersebut dapat diatasi secara dinamik dengan menempatkan sementara data partikel
pada linked-list sehingga pada struktur data node cukup menggunakan satu buah pointer. Secara umum lama waktu eksekusi dari algoritma BH ini juga bergantung pada faktor ketelitian θ yang digunakan sebagai open-test pada perhitungan gaya pada tiap partikel, dimana θ adalah fungsi dari aplikasi apa algoritma ini akan digunakan. Pada saat distribusi partikel pada sistem tidak serba sama, pada level yang sama dari tree yang terbentuk akan terdapat node kosong. Sehingga pada pembentukan tree perlu dioperasikan teknik koleksi sampah (garbage collection) tertentu.
DAFTAR PUSTAKA
1. GIBBON, P., A Hierarchical Tree Code for Quantum Molecular Dynamic Simulations, IBM Technical Report TR: 75.92.21 (1992)
2. LARS HERNQUIST, Hierarchical N-Body Methods, Computer Physics Communications 48 (1988) p. 107 – 115
3. BARNES, J. and HUT, P., A Hierarchical O(N log N) Force-Calculation Algorithm, Nature 324 (1987) p. 446 – 449
4. HANSEN, J.P., Computer Simulation of Basic Plasma Phenomena, Scottich Summer School Proceeding (1989) p. 433 – 496
5. BARNES, J. and LARS HERNQUIST, Computer Models of Colliding Galaxies, Physics Today (1993), p. 54 - 61
DISKUSI
SARIFUDDIN MADENDA 1. Dalam pencarian tree, Bapak menggunakan BFS, informasi apa yang dipakai untuk menentukan solusi untuk menentukan node yang dicari? (apakah jarak terdekat atau gaya terbesar?) 2. Bagaimana bila digunakan metode beam search?
SLAMET SANTOSA 1. Informasi acuan yang dipakai pada penjejakan tree dengan metode BFS adalah kriteria D/r ≤ θ, dimana D = lebar sektor (node) dan r = jarak dari partikel yang dihitung gayanya ke pusat node. Penjejakan dilakukan dari root, apabila kriteria tersebut belum dipenuhi, penjejakan dilakukan terus sampai node yang memenuhi atau sampai leaf. Dibandingkan dengan metode DFS (Depth First Search) metode BFS lebih efisien dan tidak mengalami penjejakan ulang (back tracking). 2. Terima kasih atas informasi tentang adanya metode beam search yang penjejakannya level demi level seperti halnya metode BFS, akan saya pelajari lebih lanjut tentang kapasitas memori yang digunakan untuk kedua metode tersebut. Kalau metode beam search adalah level demi level, dan pada prinsipnya juga tidak mengalami penjejakan ulang.
DAFTAR RIWAYAT HIDUP
1. Nama
: SLAMET SANTOSA
2. Tempat/Tanggal Lahir
: Kebumen, 23 Februari 1958
3. Instansi
: P3TM - BATAN
4. Pekerjaan / Jabatan
: Staf. Bidang Akselerator
5. Riwayat Pendidikan
: (setelah SMA sampai sekarang)
• PAT-UGM, (D3) (1984)
(D3)
• ITS, (1990)
(S1)
• Univ. of Gifu (Japan) (1995)
(S2)
6. Pengalaman Kerja
:
• 1979- 1986 : Teknisi Bidang Fisika (Instrumentasi Nuklir) • 7. Organisasi Professional
1990– sekarang : Staf Peneliti Bidang Akselerator :
• ICMI Orsat Jepang Tengah • Himpunan Fisika Indonesia Cabang Yogyakarta
HOME
KOMPUTASI DALAM SAIN DAN TEKNOLOGI NUKLIR X
