Tugas tubuh N atau cara meledakkan galaksi tanpa meninggalkan dapur

Belum lama ini saya membaca novel fiksi ilmiah "The Three-Body Problem" oleh Liu Qixin . Di dalamnya, beberapa alien memiliki masalah - mereka tidak tahu bagaimana, dengan akurasi yang cukup bagi mereka, untuk menghitung lintasan planet asal mereka. Tidak seperti kita, mereka hidup dalam sistem bintang tiga, dan "cuaca" di planet ini sangat bergantung pada lokasi relatif mereka - dari panas pembakaran hingga dingin es. Dan saya memutuskan untuk memeriksa apakah kita dapat memecahkan masalah seperti itu.

Fisika dari fenomena tersebut

Untuk memahami masalahnya, perlu berurusan dengan fisika dari fenomena tersebut. Dalam kerangka teori klasik, gaya tarik menarik dua benda ditentukan oleh hukum Newton :

$$

$$\ vec {F} (\ vec {r} _1, \ vec {r} _2) = - G m_1 m_2 \ frac {\ vec {r} _1 - \ vec {r} _2} {| \ vec {r} _1 - \ vec {r} _2 | ^ 3},$$

dimana $$$\ vec {r} _1, \ vec {r} _2$ - posisi benda di ruang angkasa, $$$m_1, m_2$ - massa tubuh, $$$G$ - Konstanta gravitasi.
Dalam sistem $$$N$ Tubuh pada masing-masing dari mereka akan dipengaruhi oleh gaya tarik dari yang lain, yang dinyatakan oleh persamaan:

$$

$$\ vec {F} _n = -G \ sum_ {k \ neq n} m_n m_k \ frac {\ vec {r} _n - \ vec {r} _k} {| \ vec {r} _n - \ vec {r } _k | ^ 3}.$$

Menggunakan hukum kedua Newton, kami menulis akselerasi untuk setiap partikel:

$$

$$\ vec {a} _n = \ vec {F} _n / m_n = -G \ sum_ {k \ neq n} m_k \ frac {\ vec {r} _n - \ vec {r} _k} {| \ vec { r} _n - \ vec {r} _k | ^ 3}.$$

Mengingat bahwa akselerasi adalah turunan kedua dari koordinat, kita memperoleh persamaan diferensial parsial orde kedua, yang harus dipecahkan untuk mendapatkan lintasan masing-masing badan:

$$

$$\ frac {\ partial ^ 2 \ vec {r} _n} {\ partial t ^ 2} = f_n = -G \ sum_ {k \ neq n} m_k \ frac {\ vec {r} _n - \ vec {r } _k} {| \ vec {r} _n - \ vec {r} _k | ^ 3}.$$

Penting untuk dicatat di sini bahwa kompleksitas fungsi komputasi $$$f_n$ sama dengan $$$O (N ^ 2)$ dan meningkat secara signifikan dengan meningkatnya jumlah badan yang berinteraksi.

Matematika

Metode pertama dan paling sederhana untuk menyelesaikan persamaan diferensial adalah metode Euler , yang dirancang untuk menyelesaikan persamaan bentuk:

$$

$$\ frac {dy} {dx} = f (x, y).$$

Setelah transisi ke wilayah diskrit, kami memperoleh:

$$

$$y_ {i} = y_ {i-1} + hf (x_ {i-1}, y_ {i-1}), \ quad i = 1,2,3, \ dots, m,$$

dimana $$$h$ Apakah langkah integrasi, dan $$$m$ - jumlah langkah integrasi. Jadi, jika kita perlu menghitung posisi tubuh pada suatu waktu $$$T$ maka kita harus melakukannya $$$m = T / h$ langkah integrasi. Di sini masalah pertama terlihat - jika $$$T$ besar, maka kita perlu mengambil sejumlah besar langkah integrasi.

Untuk menerapkan metode Euler ke masalah kita, itu harus direduksi menjadi sistem orde pertama. Untuk melakukan ini, kami memperkenalkan variabel tambahan - kecepatan partikel:

$$

$$\ frac {\ partial \ vec {v} _n} {\ partial t} = f_n, \\ \ frac {\ partial \ vec {r} _n} {\ partial t} = \ vec {v} _n.$$

Masalah kedua dalam memecahkan sistem persamaan diferensial adalah keakuratan solusi dan kontrolnya. Akurasi dapat ditingkatkan dengan dua cara: dengan mengurangi langkah integrasi dan memilih metode dengan tingkat akurasi yang lebih tinggi. Kedua metode mengarah pada peningkatan kompleksitas komputasi, tetapi dengan cara yang berbeda. Misalnya, Anda dapat menggunakan metode Runge-Kutta orde empat klasik, membutuhkan empat perhitungan fungsi $$$f_n$ di setiap langkah, tetapi memiliki urutan akurasi $$$O (h ^ 4)$ (Sebagai perbandingan, metode Euler memiliki urutan akurasi $$$O (h)$ dan membutuhkan satu perhitungan $$$f_n$ ) Keakuratan solusi juga dapat dikontrol dalam beberapa cara: membandingkan dengan solusi analitis, menyelesaikan dengan menggunakan metode yang berbeda atau dengan langkah yang berbeda dan membandingkan hasilnya, mengontrol parameter dan batasan pihak ketiga yang harus dipenuhi oleh solusi.
Juga, masing-masing metode ini memiliki kekurangan. Keputusan analitis mungkin tidak ada, atau, bahkan dalam kebanyakan kasus, sama sekali tidak ada. Misalnya, untuk tugas kita $$$N$ solusi analitik tel hanya untuk $$$N = 2$ , tetapi bahkan ini sudah cukup untuk menguji keakuratan metode. Memecahkan masalah dengan dua metode atau dengan langkah-langkah berbeda meningkatkan waktu perhitungan, tetapi pendekatan ini dapat diterapkan untuk hampir semua tugas. Tidak setiap tugas memiliki keterbatasan, tetapi untuk kita, mereka memiliki: pada setiap langkah integrasi kita dapat mengontrol implementasi undang-undang konservasi . Pendekatan ini juga meningkatkan kompleksitas perhitungan, tetapi ada banyak pilihan, menghitung jumlah momentum atau momen sudut dari semua partikel memiliki kompleksitas urutan. $$$O (N)$ , sementara menghitung energi total suatu sistem memiliki kompleksitas urutan $$$O (N ^ 2)$

Gambar 1. Contoh lintasan elips.

Pertimbangkan kasus sederhana dari satelit yang bergerak dalam orbit elips di sekitar Bumi. Ketika satelit mendekati Bumi, kecepatannya meningkat, dan ketika bergerak menjauh dari Bumi, ia menurun, karenanya, kemungkinan muncul untuk mengurangi langkah integrasi tepat waktu di bagian hijau orbit, dan meningkatkannya dalam warna merah dan biru tanpa mengubah keakuratan solusi. Mari kita coba membandingkan dengan lebih detail.

Tabel 1. Metode yang diselidiki untuk menyelesaikan persamaan diferensial

Ditandai	Memesan	Deskripsi
adams	1-5	Metode Adams-Bashfort
Euler	1	Metode Euler
rk4	4	Metode Runge-Kutta Klasik
rkck	5	Metode Kash-Karp
rkdp	5	Metode Dorman-Prince
rkdverk	6	Werner Metode 1) hal. 181
rkf	7	Felberg Metode 1) hal. 180
rkgl	6	Metode Gauss-Legendre Tersirat
rklc	4	Metode Lobatto implisit
trapeze	2	Metode trapesium

Untuk memilih metode terbaik untuk tugas kami, kami akan membandingkan beberapa metode yang diketahui. Untuk melakukan ini, kami mensimulasikan tabrakan dua sistem tubuh $$$N = 512$ dan mengukur perubahan relatif dalam total energi , momentum dan momennya pada akhir simulasi (waktu simulasi maksimum $$$T_ {maks} = 2,5$ ) Dalam hal ini, kami akan memvariasikan langkah dan parameter metode integrasi dan mengukur jumlah panggilan fungsi $$$f_n$ , masing-masing, metode-metode yang dengan jumlah panggilan yang lebih kecil akan menyebabkan lebih sedikit kerugian, kami akan anggap lebih dapat diterima.

a)	b)
Gambar 2. Perubahan relatif dalam energi a), momentum b), pada akhir simulasi sistem $$$N = 512$ tubuh dengan berbagai metode tergantung pada jumlah perhitungan fungsi $$$f_n$ presisi ganda

Dari grafik Gambar 2 dapat dilihat bahwa rasio terbaik dari jumlah perhitungan fungsi $$$f_n$ dan perubahan relatif dalam energi sistem tubuh dalam metode Adams dan Dorman-Prince orde kelima. Terlihat juga bahwa untuk semua metode dengan peningkatan jumlah perhitungan $$$f_n$ perubahan relatif dalam momentum sistem meningkat. Untuk perubahan relatif dalam energi, ini juga terlihat, tetapi hanya untuk beberapa metode yang dapat mencapai ambang batas $$$dE / E_0 <10 ^ {- 12}$ . Efek ini tidak lagi dapat dikaitkan dengan kesalahan metode, tetapi karena kesalahan perhitungan, dan peningkatan akurasi lebih lanjut hanya dimungkinkan bersama-sama dengan peningkatan akurasi perhitungan dengan floating point.

a)	b)
Gambar. 3. Perubahan relatif dalam energi a), momentum b), pada akhir simulasi sistem $$$N = 512$ tubuh dengan berbagai metode tergantung pada jumlah perhitungan fungsi $$$f_n$ dengan presisi empat kali lipat (__float128)

Gambar 3a dan 3b menunjukkan bahwa penggunaan perhitungan dengan akurasi empat kali lipat dapat mengurangi kehilangan energi relatif hingga $$$10 ^ {- 23}$ , tetapi Anda perlu memahami bahwa waktu komputasi ditingkatkan dua urutan besarnya dibandingkan dengan presisi ganda. Seperti dalam kasus perhitungan presisi ganda, rasio akurasi terbaik dan jumlah perhitungan fungsi $$$f_n$ memiliki metode urutan kelima Adams dan Dorman-Prince.

Metode Dorman-Prince dan Werner termasuk dalam kelas metode bersarang dan secara bersamaan dapat menghitung dua solusi dengan urutan akurasi tinggi dan rendah (untuk metode Dorman-Prince, pesanan 5 dan 4, dan untuk metode Werner, pesanan 6 dan 5). Jika kedua solusi ini sangat berbeda, maka kita dapat memecah langkah integrasi saat ini menjadi yang lebih kecil. Itu memungkinkan kami untuk secara dinamis mengubah langkah integrasi dan menguranginya hanya di area-area di mana diperlukan.

Mari kita bandingkan metode Dorman-Prince, Werner, dan Adams urutan kelima secara lebih rinci, selama interval simulasi yang lebih lama dari sistem kami ( $$$T_ {maks} = 300$ )

a) Perubahan relatif dalam energi, momentum dan momennya dalam proses pemodelan

b) Jumlah perhitungan fungsi $$$f_n$ pada interval simulasi $$$\ Delta T = 0,3$

Gambar 4. Ketergantungan perubahan relatif dalam parameter fisik sistem dan jumlah perhitungan fungsi $$$f_n$ dari waktu ke waktu pemodelan dengan metode variabel-langkah Dorman-Prince $$$h = 10 ^ {- 5} \ dots 10 ^ {- 9}$

a) Perubahan relatif dalam energi, momentum dan momennya dalam proses pemodelan

b) Jumlah perhitungan fungsi $$$f_n$ pada interval simulasi $$$\ Delta T = 0,3$

Gambar. 5. Ketergantungan perubahan relatif dalam parameter fisik sistem dan jumlah perhitungan fungsi $$$f_n$ versus pemodelan waktu dengan metode Werner dengan langkah variabel $$$h = 10 ^ {- 5} \ dots 10 ^ {- 9}$

Gambar. 6. Perubahan relatif dalam energi, momentum dan momennya dalam proses pemodelan dengan metode Adams-Bashfort orde lima dengan langkah $$$h = 10 ^ {- 6}$

Gambar 7. Ketergantungan jumlah perhitungan fungsi $$$f_n$ untuk metode kelima, metode Adams, Dorman-Prince dan Werner dari waktu simulasi

Dapat dilihat bahwa pada tahap awal evolusi sistem kami ( $$$T <$ 2$ ) ketiga metode tersebut menunjukkan karakteristik yang serupa, tetapi pada tahap selanjutnya peristiwa tertentu terjadi dalam sistem, akibatnya kesalahan dalam parameter utama sistem (energi total, momentum dan momentumnya) melonjak tajam. Tetapi metode Dorman-Prince dan Werner mengatasi perubahan ini jauh lebih baik karena kemampuan untuk mengurangi langkah integrasi di bagian "kompleks", sebagai akibatnya jumlah perhitungan fungsi meningkat $$$f_n$ seperti yang terlihat pada Gambar 4b dan 5b, tetapi jumlah total perhitungan $$$f_n$ metode bersarang lebih kecil dari metode Adams-Bashfort, seperti yang dapat dilihat pada Gambar 7.


Video tersebut menunjukkan hal itu hingga ke titik waktu $$$T = 25$ gerakannya relatif tenang, dan setelah tabrakan pusat "galaksi" terjadi, yang mengarah ke perubahan tajam dalam lintasan dan peningkatan kuat dalam kecepatan beberapa partikel. Selain itu, untuk menjaga keakuratan solusi, perlu untuk mengurangi langkah integrasi. Metode bersarang dapat melakukan ini secara otomatis, grafik menunjukkan bahwa di beberapa bagian evolusi sistem, langkah integrasi berkurang hampir dua urutan besarnya dengan $$$10 ^ {- 5}$ sebelumnya $$$h = 10 ^ {- 7}$ . Saat menggunakan metode Adams dan langkah seperti itu pada seluruh interval evolusi sistem, kami tidak akan menunggu solusi.

Ringkasan

Untuk solusinya, lebih baik menggunakan metode bersarang yang memungkinkan Anda mengontrol langkah integrasi secara dinamis, dan menguranginya hanya pada bagian "kompleks" dari lintasan.

Jangan mengejar metode urutan tertinggi. Bahkan ketika menggunakan tipe data 'ganda', mereka tidak mencapai kemampuan potensial mereka, menggunakan tipe data dengan akurasi yang lebih besar sangat meningkatkan waktu yang dibutuhkan untuk menyelesaikan masalah.

Implementasi CPU

Sekarang pilihan metode untuk menyelesaikan persamaan didefinisikan, mari kita coba mencari tahu perhitungan gaya interaksi untuk setiap partikel. Kami mendapatkan siklus ganda untuk semua partikel:


Gaya gravitasi untuk masing-masing tubuh dihitung secara independen, dan untuk menggunakan semua inti prosesor, cukup untuk menulis arahan OpenMP sebelum siklus pertama:


Karena Karena setiap tubuh berinteraksi dengan masing-masing, untuk mengurangi jumlah interaksi prosesor dengan RAM dan meningkatkan penggunaan cache, kami memiliki kemampuan untuk memuat bagian dari data ke dalam cache dan menggunakannya berulang kali:


Optimalisasi lebih lanjut terdiri dari mengeluarkan fungsi menghitung gaya dalam siklus utama dan menghilangkan pembagian dan perkalian dengan massa tubuh m1 [b1]. Selain fakta bahwa kami telah sedikit mengurangi perhitungan, kompiler akan dapat menerapkan instruksi vektor prosesor SSE dan AVX pada siklus yang diperluas.


Tabel 2. Ketergantungan waktu perhitungan (dalam detik) dari fungsi $$$f_n$ pada jumlah badan yang berinteraksi $$$N$ untuk berbagai implementasi CPU

$$$N$	2048	4096	8192	16384	32768
sederhana	0,0425	0,1651	0,6594	2.65	10.52
openmp	0,0078	0,0260	0,1079	0,417	1.655
openmp + blok + optimisasi	0,0037	0,0128	0,0495	0,194	0,774

Terlihat jelas bahwa versi dengan dukungan OpenMP dipercepat sebanyak enam kali, tepatnya dalam hal jumlah inti, dan versi yang dioptimalkan lebih cepat sedikit lebih dari dua kali. Jadi, dengan optimisasi, Anda seharusnya tidak hanya mengandalkan konkurensi. Menariknya, dalam perhitungan aliran tunggal, prosesor bekerja pada frekuensi 3,6 GHz, dalam versi paralel (openmp) menjatuhkan frekuensi menjadi 3,4 GHz, dan secara paralel dan dioptimalkan (openmp + blok + optimasi) turun ke 3,3 GHz, tetapi ini itu tidak menghentikannya untuk bekerja 13,6 kali lebih cepat. Terlihat juga bahwa peningkatan waktu komputasi dengan peningkatan ukuran masalah adalah kuadratik, dan lebih jauh $$$N$ membuat tugas tidak terpecahkan dalam waktu yang wajar.

Implementasi GPU

Tetapi ada keinginan untuk membuat perhitungan lebih cepat. Ada beberapa arah yang tersedia untuk akselerasi: Perhitungan GPU, perkiraan fungsi $$$f_n$ . Pertama, untuk komputasi GPU, saya memilih teknologi OpenCL. Untuk pekerjaan yang lebih nyaman, perpustakaan CLHPP digunakan. Keuntungan utama OpenCL adalah bahwa kode dapat dijalankan pada prosesor dan GPU, yang menyederhanakan penulisan dan debugging, serta memperluas daftar perangkat keras yang akan dijalankan. Alat Oclgrind membantu dalam debugging, yang pada saat runtime menunjukkan kesalahan panggilan OpenCL API dan masalah akses memori.

Opencl

Untuk memulai dengan OpenCL, Anda perlu mendapatkan daftar platform yang tersedia. Platform yang paling umum adalah AMD, Intel, dan NVidia.


Selanjutnya, setelah memilih platform, Anda harus memilih perangkat komputasi yang diwakili platform ini:


Dan pada akhir fase persiapan, perlu untuk membuat konteks dan antrian di mana memori akan dialokasikan dan perhitungan akan dilakukan. Misalnya, konteks yang menggabungkan semua perangkat komputasi dari platform yang dipilih dibuat sebagai berikut:


Untuk mengunduh kode sumber ke perangkat komputer, kode itu harus dikompilasi, kelas cl :: Program ditujukan untuk ini.


Untuk menjelaskan parameter fungsi (kernel) yang dijalankan pada perangkat komputasi, ada templat cl: make_kernel.


Lebih lanjut, semuanya sederhana: kita mendeklarasikan variabel dengan tipe kernel, meneruskan program yang dikompilasi dan nama kernel komputasi ke dalamnya, kita dapat memulai kernel hampir seperti fungsi normal.


Inti komputasi untuk OpenCL sendiri sangat mirip dengan opsi 'openmp + block + optimation' untuk CPU, hanya tidak seperti versi CPU, siklus pertama dikontrol menggunakan OpenCL (kisaran siklus ditentukan oleh variabel global_range dari kode startup kernel, dan nomor iterasi saat ini tersedia dari kernel menggunakan fungsi get_global_id (0)). Pertama, sebagian data tubuh dimuat dari memori lokal, lalu diproses. Memori lokal umum untuk semua utas dalam grup, jadi unduhan terjadi sekali untuk grup, dan diproses oleh setiap utas dalam grup dan, karena memori lokal jauh lebih cepat daripada global, perhitungannya jauh lebih cepat.

Cuda

Implementasi untuk platform NVidia CUDA sedikit lebih sederhana dari OpenCL, kami tidak perlu membuat konteks perangkat dan mengelola antrian eksekusi (setidaknya sampai kami ingin membuat implementasi multi-GPU). Seperti dalam kasus OpenCL, kita perlu mengalokasikan memori pada GPU, menyalin data kita ke dalamnya, dan kemudian kita dapat memulai inti komputasi.

Anda dapat membaca lebih lanjut tentang bekerja dengan CUDA di sini .


Tidak seperti OpenCL, CUDA tidak menentukan rentang penuh iterasi (dalam implementasi OpenCL itu adalah global_range), tetapi menetapkan ukuran kisi dan ukuran blok dalam kisi, sehingga, perhitungan jumlah tubuh saat ini sedikit berubah, jika tidak, kernel sangat mirip dengan OpenCL, dengan pengecualian nama lain fungsi sinkronisasi dan penentu untuk memori bersama. Fitur pembeda lain yang berguna dari CUDA adalah bahwa kita dapat menentukan ukuran memori bersama yang diperlukan saat kernel dimulai. Seperti dalam implementasi OpenCL, pada awal setiap blok iterasi kami menyalin bagian dari data ke memori bersama dan kemudian bekerja dengan memori ini dari semua utas blok.


Tabel 3. Ketergantungan fungsi waktu perhitungan (dalam detik) $$$f_n$ pada jumlah badan yang berinteraksi $$$N$ untuk berbagai implementasi GPU dan implementasi CPU yang lebih baik

$$$N$	4096	8192	16384	32768	65536	131072
openmp + blok + optimisasi	0,0128	0,0495	0,194	0,774	---	---
OpenCL + setengah NVidia K80	0,004	0,008	0,026	0,134	0,322	1.18
CUDA + setengah NVidia K80	0,004	0,008	0,0245	0,115	0,291	1.13

Secara umum, hasilnya mengecewakan, hanya 5-6 kali lebih cepat daripada implementasi CPU. Bahkan jika kita melakukan perhitungan pada seluruh K80, kita akan mendapatkan akselerasi hingga 12 kali, tetapi sejak itu Karena kompleksitas tugasnya kuadratik, maka dalam waktu yang masuk akal kita akan dapat memproses bukan 32.768 badan yang berinteraksi, tetapi 131072, yang hanya 4 kali lebih banyak.

Perkiraan fungsi $$$f_n$

Jika Anda melihat dekat pada fungsi, yang diatur oleh kekuatan menarik dari dua tubuh, jelas bahwa itu berkurang secara kuadrat dengan jarak. Oleh karena itu, kita dapat secara akurat menghitung kekuatan interaksi antara benda yang dekat, dan kira-kira antara yang jauh. Satu pendekatan terkenal
adalah algoritma treecode yang diusulkan oleh D. Barnes dan P. Hat. Pohon octree dibangun di ruang simulasi, berisi koordinat daun dan massa badan simulasi. Simpul induk mengandung pusat massa, massa total simpul anak dan jari-jari bola yang dijelaskan di sekitar tubuh simpul anak. Akar pohon berisi pusat massa semua benda, massa totalnya, dan jari-jari bola yang dijelaskan di sekelilingnya. Saat menghitung gaya interaksi, jarak ke akar pohon pertama kali dipertimbangkan jika rasio jarak ke simpul ke jari-jarinya lebih besar dari beberapa konstanta $$$\ lambda_ {crit}$ , maka root dianggap satu tubuh dengan koordinat sama dengan koordinat pusat massa tubuh yang termasuk di dalamnya, dan massa sama dengan jumlah massa simpul anak perempuan, jika rasionya kurang dari atau sama dengan $$$\ lambda_ {crit}$ , lalu prosedur ini diulang secara rekursif untuk setiap simpul anak. Jika daun pohon tercapai, maka kekuatan interaksi dipertimbangkan dengan cara biasa. Jadi, jika satu benda jauh dari kelompok kompak dari badan lain, kelompok ini disajikan sebagai satu tubuh, dan gaya interaksi dihitung tidak sampai masing-masing tubuh, tetapi hanya hingga satu tubuh. Karena ini, kompleksitas algoritma menurun dengan $$$O (N ^ 2)$ sebelumnya $$$O (N \ log (N))$ dengan mengorbankan beberapa kehilangan akurasi.

Dalam pendekatan saya, saya memutuskan untuk tidak menggunakan pohon octree , tetapi pohon kd, karena lebih mudah digunakan dan memiliki overhead penyimpanan yang lebih rendah dibandingkan dengan octree.

Mari kita kembali ke implementasi pada CPU. Node kd-tree dapat direpresentasikan dalam bentuk kelas yang berisi pointer ke turunan kiri dan kanan dan informasi tentang koordinat dan massa:


Dengan metode ini menyimpan pohon, kami memiliki dua opsi yang mungkin untuk melintasi pohon: baik menggunakan rekursi eksplisit, atau gunakan tumpukan sendiri. Saya memilih opsi kedua.


Seperti dalam kasus implementasi CPU "tepat", fungsi perhitungan gaya dipanggil untuk setiap tubuh. Siklus di semua badan dapat dengan mudah diparalelkan menggunakan arahan OpenMP.

Tetapi iterasi loop tetangga dalam kasus ini akan merujuk ke bagian pohon yang sama sekali berbeda, yang tidak memungkinkan penggunaan cache prosesor secara efisien. Untuk mengatasi masalah ini, semua benda harus dilintasi bukan dalam urutan asli, tetapi dalam urutan di mana mayat berada di daun pohon kd, maka iterasi yang berdekatan akan terjadi untuk benda yang dekat di ruang angkasa, dan akan melintasi pohon di sepanjang jalur yang hampir sama.


Implementasi ini memiliki masalah lain - tidak ada cara universal untuk memaralelkan traversal pohon semacam itu. Tetapi kita dapat sepenuhnya mengubah cara pohon disimpan dalam memori - kita dapat menyimpan semua node dalam satu array linier dan sepenuhnya meninggalkan penyimpanan pointer ke keturunan, dengan analogi dengan pembangunan tumpukan biner . Saat memulai penomoran node dengan $$$1$ jika simpul ada di nomor sel $$$i$ , lalu anak kirinya ada di sel $$$2i$ , anak yang tepat di sel $$$2i + 1$ , dan induk di sel $$$i / 2$ . Node kanan sesuai dengan kiri dengan nomor $$$i$ akan memiliki nomor $$$i + 1$ . Array itu sendiri akan memiliki panjang $$$2N$ , dan semua simpul daun akan ditemukan di bagian terakhir $$$N$ sel, apalagi, node yang berjarak dekat akan berlokasi di sel dekat array. Fungsi traversal pohon akan berubah sedikit, tetapi kerangka kerjanya tetap sama:


Tapi ini tidak semua kemungkinan bahwa penyimpanan node dalam array terbuka untuk kita - kita bisa menolak stack saat merangkak. Untuk melakukan ini, dalam cabang kode di mana kita pergi ke anak-anak dari simpul, kita menambahkan fungsi untuk menghitung simpul berikutnya ( $$$next_ {up}$ ), dan di cabang tempat kami menghitung gaya interaksi, kami menambahkan perhitungan simpul berikutnya dengan subtree saat ini dilewati ( $$$skip_ {down}$ )

Untuk melewati subtree saat ini, kita perlu turun ke root (arah $$$D$ ), ketika kita berada di keturunan kanan, segera setelah kita mencapai kiri, kita pergi ke subtree kanan yang sesuai (arah $$$R$ ), jika kita sampai ke root, maka traversal pohon selesai.



Gambar 8. Melewati subtree $$$skip_ {down}$ .

Untuk pergi ke simpul berikutnya, jika mungkin, pergi ke anak kiri (arah $$$U$ ), dan jika tidak ada turunan, maka pergi ke simpul berikutnya 'dari bawah' menggunakan fungsi $$$skip_ {down}$ .



Gambar 9. Transisi ke node berikutnya $$$next_ {up}$ .

Tampaknya kami mengganti rekursi dengan loop $$$sementara$ dalam fungsi $$$skip_ {down}$ , dan penggantian seperti itu tidak memberikan apa-apa, tetapi mari kita lihat bagaimana menentukan apakah suatu simpul dengan angka $$$i$ keturunan benar. Ini dapat dilakukan hanya dengan memeriksa nomor ganjilnya (anak yang tepat ada di dalam sel $$$2i + 1$ ), untuk ini cukup untuk menghitung $$$i \ & 1$ . Yaitu kami membagi nomornya $$$i$ dua jika bit paling signifikan diatur ke satu. Tapi ini bisa dilakukan tanpa loop, di banyak prosesor ada instruksi set pertama yang mengembalikan posisi bit set pertama, dengan menggunakannya kita mengubah loop menjadi empat instruksi prosesor.


Setelah itu kita bisa mengecualikan stack dari fungsi tree traversal dan menggantinya dengan pair $$$skip_ {down} + next_ {up}$ , setelah itu fungsinya terlihat lebih sederhana:


Secara total, kami mendapat enam kombinasi traversal pohon dan perhitungan gaya. Bandingkan pendekatan ini dalam hal waktu dan kualitas perhitungan. Kami mengambil sebagai ukuran kualitas perubahan relatif dalam total energi sistem setelah 100 iterasi. Sebagai sistem model, kami mengambil dua "galaksi" yang berinteraksi yang terdiri dari$$$16384$tubuh masing-masing.

Tabel 4. Kombinasi metode traversal pohon dan perhitungan gaya

Tipe perhitungan traversal / gaya pohon	Pohon dengan tumpukan	'Tumpukan' dengan tumpukan	'Tumpukan' tanpa tumpukan
Iterasi berdasarkan nomor tubuh	cycle + tree	cycle + heap	cycle + heapstackless
Bypass daun	nestedtree + tree	nestedtree + heap	nestedtree + heapstackless

a) Fungsi perhitungan ketergantungan waktu $$$f_n$ dari rasio jarak kritis ke simpul pohon ke jari-jarinya ( $$$\lambda_{crit}$ )	b) Ketergantungan perubahan relatif energi dalam sistem pada rasio jarak kritis ke simpul pohon terhadap jari-jarinya ( $$$\lambda_{crit}$ )
Gambar 10. Hasil perhitungan fungsi $$$f_n$ dalam berbagai cara pada CPU (jumlah tubuh $$$N=32768$ )

Dapat dilihat bahwa implementasi 'nestedtree + tree' tidak ada harapan dalam kecepatan, karena tidak memiliki konkurensi. Implementasi dengan lokasi node pohon dalam array dan pengindeksan seperti dalam tumpukan biner berada di depan. Perubahan relatif dalam energi dapat diabaikan untuk semua varian dengan$$$\lambda_{crit} > 1$ .Juga dalam gambar. 10a menunjukkan bahwa untuk ($$$\lambda_{crit} < 10$ ) semua varian perhitungan fungsi $$$f_n$ secara signifikan menyalip dalam kecepatan versi yang paling optimal dari perhitungan yang tepat ('openmp + block + optimization'), dengan peningkatan lebih lanjut $$$\lambda_{crit}$ implementasi dengan pohon kehilangan versi pastinya.

Traversal pohon GPU

Saya mencoba memintas pohon pada GPU menggunakan teknologi OpenCL dan CUDA. Opsi menyimpan node dalam bentuk pohon segera dibuang, dan hanya opsi yang tersisa dengan menyimpan pohon dalam array dengan pengindeksan seperti dalam tumpukan biner. Secara umum, implementasi core komputasi tidak jauh berbeda dengan versi CPU.


Pada versi pertama, tree traversal dimulai dengan urutan penomoran tubuh pada array asli, sehingga utas tetangga melewati bagian-bagian pohon yang sama sekali berbeda, yang secara negatif mempengaruhi kinerja cache memori GPU. Oleh karena itu, dalam perwujudan kedua, sebuah traversal diterapkan, mulai dari atas pohon, dalam hal ini, benang tetangga mulai melintasi pohon di sepanjang jalan yang sama, karena puncak pohon tetangga dekat dan di luar angkasa. Penting juga bahwa kita memilih penomoran dalam array node pohon bukan dari awal, tetapi dari satu, dalam hal ini daun pohon disimpan di paruh kedua array, dan dengan jumlah tubuh yang sama dengan kekuatan dua, kita akan memiliki akses yang sama ke memori dengan indeks tn1.


Saat melintasi urutan penomoran simpul pohon, kami mendapat peningkatan kinerja. Tetapi opsi ini juga dapat ditingkatkan. Memori global tempat node pohon saat ini dioptimalkan untuk akses bersama , mis. utas satu grup harus membaca kata-kata yang terletak di satu blok memori. Dalam kasus kami, traversal pohon dimulai di sepanjang jalur yang sama, dan kami meminta data yang sama dengan semua utas grup, tetapi semakin jauh kami masuk lebih dalam ke pohon, jalur benang tetangga semakin berbeda, dan kami harus meminta data yang berbeda, yang mengurangi kinerja subsistem beberapa kali memori. Tetapi node dari masing-masing subtree milik tingkat yang sama terletak di sel-sel memori yang relatif dekat. Yaituketika melintasi sisa pohon, utas tetangga dari inti komputasi tidak mengakses simpul pohon yang sama, tetapi berada dalam memori. Untuk mengoptimalkan akses memori ini, memori tekstur dapat digunakan. Tapi ada satu halangan. Saat ini, tekstur tidak mendukung kerja dengan data presisi ganda (kami ingin menghitung secara akurat). Tetapi dalam CUDA ada fungsi __hiloint2double yang mengumpulkan angka presisi ganda dari dua bilangan bulat.


Dalam hal ini, dua implementasi dibuat, dalam satu setiap elemen pohon (x, y, z, tree_crit_r2) diminta secara independen, dan dalam implementasi kedua permintaan ini digabungkan. Permintaan massa node terjadi jauh lebih jarang, hanya jika kondisi r2> tree_crit_r2 [curr] terpenuhi , jadi tidak masuk akal untuk menggabungkan permintaan ini dengan yang lain. Fitur lain yang berguna dari kerangka kerja CUDA adalah kemampuan untuk mengontrol rasio ukuran cache L1 dan ukuran memori bersama ( cudaFuncSetCacheConfig ). Dalam hal traversal pohon, kami tidak menggunakan memori bersama, sehingga kami dapat meningkatkan cache L1 hingga merusaknya.


Analisis program dalam profiler nvprof menunjukkan bahwa bahkan ketika menggunakan memori tekstur untuk menyimpan pohon, masih ada beban yang sangat tinggi pada memori global.

Memang, dalam CUDA, semua memori kernel yang ditujukan ke alamat 'dihitung' disimpan dalam memori global, dan karenanya, tumpukan yang diperlukan untuk melintasi pohon terletak di memori global dan menggerogoti bagian signifikan dari bandwidth chip memori, karena tumpukan tersebut memiliki masing-masing menjalankan thread, dan ada banyak utas.

Tapi, untungnya, kita sudah tahu cara melintasi pohon tanpa menggunakan tumpukan. Melengkapi inti komputasi sebelumnya dengan fungsi komputasi simpul pohon berikutnya, kami mendapatkan kernel baru, lebih dari itu, lebih ringkas.


Kinerja core yang berjalan pada GPU sangat tergantung pada ukuran blok yang digunakan untuk membagi tugas. Tergantung pada ukuran ini berapa banyak register, memori lokal dan sumber daya lainnya akan tersedia untuk setiap utas komputasi. Anda juga perlu mengingat bahwa sambil menunggu akses memori di satu utas, utas lain dapat melakukan perhitungan pada prosesor shader, dengan demikian, dengan jumlah yang cukup dari utas yang dijalankan secara bersamaan pada satu prosesor, waktu akses memori akan disembunyikan di balik perhitungan. Karena itu, sebelum membandingkan kinerja core kami, kami perlu menghitung ukuran blok optimal untuk masing-masing core. Mari kita bandingkan setengah yang tersedia untuk kita dari NVidia K80.

Tabel 5. Ketergantungan fungsi waktu perhitungan (dalam detik) $$$f_n$ GPU $$$N=131072$ $$$\lambda_{crit}=10$

/	8	16	32	64	128	256	512	1024
opencl+dense	5.77	2.84	1.46	1.13	1.15	1.14	1.14	1.13
cuda+dense	5.44	2.55	1.27	0.96	0.97	0.99	0.99	-
opencl+heap+cycle	5.88	5.65	5.74	5.96	5.37	5.38	5.35	5.38
opencl+heap+nested	4.54	3.68	3.98	5.25	5.46	5.41	5.48	5.31
cuda+heap+nested	3.62	2.81	2.68	4.26	4.84	4.75	4.8	4.67
cuda+heap+nested+tex	2.6	1.51	0.912	0.7	1.85	1.75	1.69	1.61
cuda+heap+nested+tex+stackless	2.3	1.29	0.773	0.5	0.51	0.52	0.52	0.52

6. ( ) $$$f_n$ GPU $$$N=1M$ $$$\lambda_{crit}=4$

/	8	16	32	64	128	256	512	1024
opencl+dense	366	179	89.9	69.3	70.3	69.1	68.9	68.0
cuda+dense	346	162	79.6	60.8	60.8	60.7	59.6	-
opencl+heap+cycle	16.2	18.2	20.1	21.2	21.2	21.3	21.2	21.1
opencl+heap+nested	10.5	7.63	6.38	8.23	9.95	9.89	9.65	9.66
cuda+heap+nested	8.67	6.44	5.39	5.93	8.65	8.61	8.41	8.27
cuda+heap+nested+tex	6.38	3.57	2.13	1.44	3.56	3.46	3.30	3.29
cuda+heap+nested+tex+stackless	5.78	3.19	1.83	1.21	1.11	1.10	1.11	1.13

Situasi yang sulit, tetapi, tidak seperti versi CPU dari traversal pohon, jelas bahwa setiap langkah optimasi membawa hasil nyata. Implementasi 'opencl + heap + cycle' hampir 6 kali lebih lambat dari solusi yang tepat dengan perhitungan fungsi penuh$$$f_n$ .Implementasi 'opencl + heap + nested', di mana traversal pohon dimulai dari node tetangga, 1,4 kali lebih cepat dari yang sebelumnya, karena Memori cache digunakan lebih baik. Dalam implementasi 'cuda + heap + nested', L1 cache ditingkatkan hingga merusak memori bersama, yang meningkatkan kecepatannya sebesar 1,4 kali, meskipun ada kemungkinan bahwa dalam implementasi cuda kernel kompilasi lebih optimal dikompilasi (dalam versi 'opencl + padat' dan 'cuda + core padat identik, dan kinerja versi cuda ~ 1,2 kali lebih tinggi). Peningkatan kecepatan komputasi yang paling signifikan (3,8 kali) dicapai ketika pohon terletak di memori tekstur dan pertanyaan ke elemen simpul pohon digabungkan. Implementasi dengan tree traversal tanpa menggunakan tumpukan 'cuda + heap + nested + tex + stackless' juga 1,4 kali lebih cepat.di dalamnya, seluruh bandwidth bus memori hanya digunakan untuk mengakses data tentang simpul pohon dan tidak dihabiskan di stack. Demikianlah dengan$$$\lambda_{crit}=10$ Dimungkinkan untuk mencapai akselerasi dua kali dibandingkan dengan perhitungan penuh fungsi $$$f_n$ . Tapi $$$\lambda_{crit}=10$ nilai parameter yang terlalu besar, pada grafik perubahan relatif energi dalam sistem dari rasio jarak kritis ke simpul pohon hingga radiusnya untuk implementasi CPU, dapat dilihat bahwa nilai yang lebih rendah dapat digunakan $$$\lambda_{crit}$tanpa kehilangan akurasi. Mari kita coba bervariasi$$$\lambda_{crit}$ dengan ukuran blok optimal yang kami tentukan pada langkah sebelumnya.

a) $$$N=128K$	b) $$$N=1M$
Gambar 11. Ketergantungan waktu perhitungan fungsi $$$f_n$ dari rasio jarak kritis ke simpul pohon ke jari-jarinya ( $$$\lambda_{crit}$ ) untuk berbagai implementasi GPU tree walk

Dapat dilihat bahwa untuk kecil $$$\lambda_{crit}$ semua metode perhitungan fungsi $$$f_n$pergi untuk menutup nilai, ditentukan oleh waktu membangun kd-tree dan menyiapkan data untuk GPU. Selain itu, waktu membangun pohon memberikan kontribusi yang signifikan terhadap total waktu$$$\lambda_{crit}\leq 4$, maka saat ini bisa diabaikan. Sangat menarik untuk dicatat bahwa kapan$$$N=128K$ kinerja meningkat lagi setelah mencapai $$$\lambda_{crit}=1024$kemungkinan besar ini disebabkan oleh kenyataan bahwa semua utas GPU mengelilingi simpul pohon yang sama pada saat yang sama, yang meningkatkan penggunaan cache L1 dan sepenuhnya menghilangkan 'cabang divergensi' . Terlihat juga bahwa kinerja terbaik untuk implementasi tanpa stack (cuda + heap + nested + tex + stackless), ini mengungguli versi dengan stack sekitar$$$1.4-1.5$kali. Implementasi lain beberapa kali lebih lambat. Untuk mengkonsolidasikan hasilnya, kami akan menghitung waktu pada GPU dengan arsitektur yang lebih baru.


Saat menggunakan GeForce GTX 1080 Ti untuk perhitungan, perbedaan antara jalan-jalan pohon dengan tumpukan dan tanpa tumpukan mencapai dua kali, asalkan kita mengabaikan waktu yang diperlukan untuk membangun pohon. Fakta ini mendorong kami untuk memastikan bahwa di masa depan porting ke GPU dan membangun pohon. Perbandingan tabel 5-7 menunjukkan bahwa tidak ada ukuran blok optimal tunggal untuk jumlah tubuh yang berbeda, dan bahkan lebih lagi untuk arsitektur GPU yang berbeda, apalagi, perbedaan waktu perhitungan dapat mencapai beberapa kali, bahkan jika Anda tidak memperhitungkan nilai batas. Jadi, sebelum perhitungan panjang, masuk akal untuk menentukan ukuran blok optimal untuk setiap konfigurasi.

Hal utama yang kami capai adalah kemampuan untuk memodelkan interaksi berpasangan dari sedikit lebih dari satu juta tubuh ($$$2^{20}$) untuk waktu yang masuk akal, bukan GPU terbaru. Pada GPU yang lebih baru (Tesla V100), tampaknya, akan mungkin untuk memproses sekitar dua juta badan yang berinteraksi dalam waktu yang wajar, karena kinerja mereka sekitar empat kali lebih tinggi daripada setengah dari Tesla K80. Meskipun jumlah ini juga tak tertandingi dengan jumlah bintang bahkan di galaksi kerdil, seperti Awan Magellan Kecil , tetapi, menurut saya, sangat mengesankan.

Kesimpulan

Penggunaan metode embedded untuk menyelesaikan persamaan diferensial memungkinkan untuk menyelesaikan masalah serupa dengan tingkat akurasi tinggi dengan biaya waktu yang relatif kecil, dan penggunaan algoritma aproksimasi untuk fungsi menghitung gaya tarik berpasangan memungkinkan kita untuk memproses volume kolosal dari benda yang saling berinteraksi. Dengan demikian, tidak seperti alien dari "Tugas Tiga Tubuh", kami dapat menyelesaikan masalah$$$N$ tubuh, dan untuk sejumlah kecil tubuh dan untuk seluruh gugus bintang, meskipun dengan mengorbankan sejumlah akurasi.

Visualisasi

Bagi mereka yang telah membaca sampai akhir, saya akan memberikan beberapa video lagi dengan visualisasi evolusi sistem tubuh.

Simulasi tabrakan dua galaksi. Jumlah mayat adalah 60 ribu.


Memodelkan evolusi sebuah galaksi yang terdiri dari sejuta bintang. Di tengah adalah tubuh dengan massa 99% dari total. Partikel tunggal hampir tidak bisa dibedakan. Sudah lebih seperti gelombang dalam setetes cairan. Berwarna sesuai dengan modul kecepatan. Kecepatan rendah - biru, sedang - hijau, tinggi - merah. Terlihat bahwa di bagian tengah kecepatannya lebih tinggi, dan secara bertahap berkurang ke tepinya, dan yang terendah di bidang ekuator.


Simulasi yang lebih 'dinamis' dari evolusi galaksi jutaan bintang. Di tengah adalah tubuh dengan massa 10% dari total. Tubuh pusat mempengaruhi sisanya secara signifikan lebih sedikit. Pertama, "bintang-bintang" terbang terpisah, kemudian berkumpul kembali menjadi beberapa kelompok, dan pada akhirnya kembali membentuk satu kelompok besar.

Dalam proses pemodelan, sekitar 5% dari "bintang" meninggalkan wilayah awal tidak dapat ditarik kembali.


Pada detik ke-10, sangat mirip dengan galaksi roda satu yang ada .

Kode proyek dapat ditemukan di github .