A tarefa de N corpos ou como explodir uma galáxia sem sair da cozinha

Não faz muito tempo, li o romance de ficção científica "O Problema dos Três Corpos", de Liu Qixin . Nele, alguns alienígenas tinham um problema - eles não sabiam como, com precisão suficiente, calcular a trajetória de seu planeta natal. Ao contrário de nós, eles viviam em um sistema de três estrelas, e o “clima” no planeta dependia fortemente de sua localização relativa - do calor incinerador ao frio gélido. E decidi verificar se podemos resolver esses problemas.

Física do fenômeno

Para entender o problema, é necessário lidar com a física do fenômeno. No quadro da teoria clássica, a força atrativa de dois corpos é determinada pela lei de 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},$$

onde $$$\ vec {r} _1, \ vec {r} _2$ - a posição dos corpos no espaço, $$$m_1, m_2$ - massas de corpos, $$$G$ - constante gravitacional.
No sistema de $$$N$ Os corpos de cada um deles serão afetados pela força de atração do resto, que é expressa pela equação:

$$

$$\ 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}.$$

Usando a segunda lei de Newton, escrevemos a aceleração para cada partícula:

$$

$$\ 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}.$$

Lembrando que a aceleração é a segunda derivada da coordenada, obtemos uma equação diferencial parcial de segunda ordem, que deve ser resolvida para obter a trajetória de cada corpo:

$$

$$\ frac {\ parcial ^ 2 \ v {{}} n \ {parcial 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}.$$

É importante notar aqui que a complexidade de calcular uma função $$$f_n$ é igual a $$$O (N ^ 2)$ e aumenta significativamente com o aumento do número de corpos em interação.

Matemática

O primeiro e mais simples método para resolver equações diferenciais é o método de Euler , projetado para resolver equações da forma:

$$

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

Após a transição para a região discreta, obtemos:

$$

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

onde $$$h$ A etapa de integração e $$$m$ - o número de etapas de integração. Assim, se precisarmos calcular a posição dos corpos de cada vez $$$T$ então devemos fazer $$$m = T / h$ etapas de integração. Aqui o primeiro problema é visível - se $$$T$ grande, precisamos dar um grande número de etapas de integração.

Para aplicar o método Euler ao nosso problema, ele deve ser reduzido a um sistema de primeira ordem. Para fazer isso, introduzimos uma variável adicional - velocidade das partículas:

$$

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

O segundo problema na solução de sistemas de equações diferenciais é a precisão da solução e seu controle. A precisão pode ser aprimorada de duas maneiras: diminuindo a etapa de integração e escolhendo um método com uma ordem de precisão mais alta. Ambos os métodos levam ao aumento da complexidade computacional, mas de maneiras diferentes. Por exemplo, você pode usar o método Runge-Kutta clássico de quarta ordem ; ele requer quatro cálculos de função $$$f_n$ a cada passo, mas tem uma ordem de precisão $$$O (h ^ 4)$ (para comparação, o método de Euler tem uma ordem de precisão $$$O (h)$ e requer um cálculo $$$f_n$ ) A precisão da solução também pode ser controlada de várias maneiras: compare com a solução analítica, resolva usando métodos diferentes ou com etapas diferentes e compare os resultados, controle parâmetros e limitações de terceiros que a solução deve cumprir.
Além disso, cada um desses métodos tem suas desvantagens. As decisões analíticas podem estar ausentes ou, mesmo na maioria dos casos, completamente ausentes. Por exemplo, para nossa tarefa $$$N$ solução analítica tel é apenas para $$$N = 2$ , mas mesmo isso é suficiente para testar a precisão dos métodos. A resolução de um problema por dois métodos ou com etapas diferentes aumenta o tempo de computação, mas essa abordagem pode ser aplicada a quase qualquer tarefa. Nem todas as tarefas têm limitações, mas as nossas têm: em cada etapa da integração, podemos controlar a implementação das leis de conservação . Essa abordagem também aumenta a complexidade do cálculo, mas há muito por onde escolher; o cálculo da soma do momento ou do momento angular de todas as partículas tem complexidade da ordem $$$O (N)$ , enquanto o cálculo da energia total de um sistema tem complexidade de pedidos $$$O (N ^ 2)$

Fig 1. Um exemplo de uma trajetória elíptica.

Considere o caso simples de um satélite se movendo em uma órbita elíptica ao redor da Terra. À medida que o satélite se aproxima da Terra, sua velocidade aumenta e, ao se afastar da Terra, diminui, consequentemente, a possibilidade de diminuir o passo de integração no tempo na parte verde da órbita e aumentá-lo em vermelho e azul sem alterar a precisão da solução. Vamos tentar comparar com mais detalhes.

Tabela 1. Métodos investigados para resolver equações diferenciais

Marcado	Encomendar	Descrição do produto
adams	1-5	Método Adams-Bashfort
euler	1	Método Euler
rk4	4	Método Runge-Kutta clássico
rkck	5	Método Kash-Karp
rkdp	5	Método Dorman-Prince
rkdverk	6	Método de Werner 1) p. 181
rkf	7	Método Felberg 1) p. 180
rkgl	6	Método implícito de Gauss-Legendre
rklc	4	Método implícito de Lobatto
trapézio	2	Método trapézio

Para selecionar o melhor método para nossa tarefa, compararemos vários métodos conhecidos. Para isso, simulamos a colisão de dois sistemas de corpos $$$N = 512$ e medir a mudança relativa na energia total , momento e seu momento no final da simulação (tempo máximo de simulação $$$T_ {max} = 2,5$ ) Nesse caso, variaremos a etapa e os parâmetros dos métodos de integração e mediremos o número de chamadas de função $$$f_n$ , respectivamente, aqueles métodos que com um número menor de chamadas levarão a menos perdas, consideraremos mais aceitáveis.

a)	b)
Figura 2. Mudança relativa na energia a), momento b), no final da simulação do sistema $$$N = 512$ corpos por vários métodos, dependendo do número de cálculos de funções $$$f_n$ precisão dupla

A partir dos gráficos da Figura 2, é possível observar que a melhor proporção da quantidade de cálculo da função $$$f_n$ e a mudança relativa na energia do sistema dos corpos nos métodos de Adams e Dorman-Prince de quinta ordem. Também é visto que para todos os métodos com um aumento no número de cálculos $$$f_n$ a mudança relativa no momento do sistema aumenta. Para uma mudança relativa na energia, isso também é perceptível, mas apenas para alguns métodos que podem atingir o limite $$$dE / E_0 <10 ^ {- 12}$ . Esse efeito não pode mais ser atribuído ao erro dos métodos, mas ao erro dos cálculos, e um aumento adicional na precisão é possível apenas junto com um aumento na precisão dos cálculos com ponto flutuante.

a)	b)
Fig. 3. Mudança relativa na energia a), momento b), no final da simulação do sistema $$$N = 512$ corpos por vários métodos, dependendo do número de cálculos de funções $$$f_n$ com precisão quádrupla (__float128)

As Figuras 3a e 3b mostram que o uso de cálculos com precisão quádrupla pode reduzir as perdas de energia relativas em até $$$10 ^ {- 23}$ , mas você precisa entender que o tempo de computação é aumentado em duas ordens de magnitude em comparação com a precisão dupla. Como no caso de cálculos de dupla precisão, a melhor razão de precisão e o número de cálculos de funções $$$f_n$ possuem métodos de Adams de quinta ordem e Dorman-Prince.

Os métodos Dorman-Prince e Werner pertencem à classe de métodos aninhados e podem calcular simultaneamente duas soluções com alta e baixa ordem de precisão (para o método Dorman-Prince, ordens 5 e 4, e para o método Werner, ordens 6 e 5). Se essas duas soluções forem muito diferentes, poderemos dividir a etapa atual de integração em outras menores. Isso nos permite alterar dinamicamente a etapa de integração e reduzi-la apenas nas áreas em que é necessária.

Vamos comparar os métodos de quinta ordem Dorman-Prince, Werner e Adams com mais detalhes, durante um intervalo de simulação mais longo do nosso sistema ( $$$T_ {max} = 300$ )

a) A mudança relativa de energia, momento e seu momento no processo de modelagem

b) Número de cálculos de funções $$$f_n$ no intervalo de simulação $$$\ Delta T = 0,3$

Fig. 4. Dependências da mudança relativa nos parâmetros físicos do sistema e o número de cálculos de funções $$$f_n$ Modelagem periódica com o método de etapas variáveis ​​de Dorman-Prince $$$h = 10 ^ {- 5} \ pontos 10 ^ {- 9}$

a) A mudança relativa de energia, momento e seu momento no processo de modelagem

b) Número de cálculos de funções $$$f_n$ no intervalo de simulação $$$\ Delta T = 0,3$

Fig. 5. Dependências da mudança relativa nos parâmetros físicos do sistema e o número de cálculos de funções $$$f_n$ modelagem versus tempo pelo método Werner com uma etapa variável $$$h = 10 ^ {- 5} \ pontos 10 ^ {- 9}$

Fig. 6. Mudança relativa de energia, momento e seu momento no processo de modelagem pelo método de Adams-Bashfort de quinta ordem com um passo $$$h = 10 ^ {- 6}$

Fig 7. Dependências do número de cálculos de funções $$$f_n$ para métodos Adams, Dorman-Prince e Werner de quinta ordem a partir do tempo de simulação

Pode-se ver que no estágio inicial da evolução do nosso sistema ( $$$T <25$ ) os três métodos mostram características semelhantes, mas em estágios posteriores alguns eventos ocorrem no sistema, como resultado dos erros nos principais parâmetros do sistema (energia total, momento e seu momento) saltam acentuadamente. Mas os métodos Dorman-Prince e Werner lidam com essas mudanças muito melhor devido à capacidade de reduzir a etapa de integração em seções "complexas", como resultado do qual o número de cálculos de funções aumenta $$$f_n$ como visto nas figuras 4b e 5b, mas o número total de cálculos $$$f_n$ métodos aninhados são menores que o método Adams-Bashfort, como pode ser visto na Figura 7.


O vídeo demonstra que até o momento $$$T = 25$ o movimento é relativamente calmo e após uma colisão dos centros das “galáxias”, o que leva a uma mudança acentuada nas trajetórias e um forte aumento nas velocidades de algumas partículas. Além disso, para manter a precisão da solução, é necessário reduzir a etapa de integração. Os métodos aninhados podem fazer isso automaticamente; os gráficos mostram que em algumas partes da evolução do sistema, a etapa de integração foi reduzida em quase duas ordens de magnitude com $$$10 ^ {- 5}$ antes $$$h = 10 ^ {- 7}$ . Ao usar o método Adams e essa etapa em todo o intervalo da evolução do sistema, não teríamos esperado por uma solução.

Sumário

Para a solução, é melhor usar métodos aninhados que permitem controlar dinamicamente a etapa de integração e reduzi-la apenas em seções "complexas" da trajetória.

Não siga os métodos de ordem mais alta. Mesmo ao usar o tipo de dados 'double', eles não atingem seus recursos potenciais, usar os tipos de dados com maior precisão aumenta muito o tempo necessário para resolver o problema.

Implementação de CPU

Agora que a escolha de um método para resolver equações está definida, vamos tentar descobrir o cálculo da força de interação para cada partícula. Temos um ciclo duplo para todas as partículas:


As forças gravitacionais de cada corpo são calculadas independentemente e, para usar todos os núcleos do processador, basta escrever a diretiva OpenMP antes do primeiro ciclo:


Porque Como cada corpo interage com cada um deles, para reduzir o número de interações do processador com a RAM e melhorar o uso do cache, temos a capacidade de carregar parte dos dados no cache e usá-lo repetidamente:


Uma otimização adicional consiste em retirar o conteúdo da função de calcular a força no ciclo principal e eliminar a divisão e multiplicação pela massa corporal m1 [b1]. Além do fato de termos reduzido um pouco os cálculos, o compilador poderá aplicar instruções vetoriais do processador SSE e AVX em um ciclo tão expandido.


Tabela 2. Dependência do tempo de cálculo (em segundos) da função $$$f_n$ no número de corpos interagindo $$$N$ para várias implementações de CPU

$$$N$	2048	4096	8192	16384	32768
simples	0,0425	0,1651	0,6594	2,65	10,52
openmp	0,0078	0,0260	0,1079	0,417	1.655
openmp + block + otimização	0,0037	0,0128	0,0495	0,194	0,774

É claramente visto que a versão com suporte ao OpenMP é acelerada em seis vezes, exatamente em termos de número de núcleos, e a versão otimizada é mais rápida em pouco mais de duas vezes. Portanto, com a otimização, você não deve confiar apenas no paralelismo. Curiosamente, nos cálculos de fluxo único, o processador trabalhava com uma frequência de 3,6 GHz, na versão paralela (openmp) reduzia a frequência para 3,4 GHz e, no paralelo e otimizado (openmp + bloco + otimização), reduzia para 3,3 GHz, mas isso isso não a impediu de trabalhar 13,6 vezes mais rápido. Também é visto que um aumento no tempo de computação com um aumento no tamanho do problema é quadrático e um aumento adicional $$$N$ torna a tarefa insolúvel em um tempo razoável.

Implementação de GPU

Mas há um desejo de fazer cálculos ainda mais rápidos. Existem várias instruções disponíveis para aceleração: computação em GPU, aproximação de funções $$$f_n$ . Primeiro, para a computação da GPU, escolhi a tecnologia OpenCL. Para um trabalho mais conveniente, a biblioteca CLHPP foi usada. A principal vantagem do OpenCL é que o código pode ser executado no processador e na GPU, o que simplifica a gravação e a depuração, além de expandir a lista de hardware a ser executada. A ferramenta Oclgrind ajuda na depuração, que em tempo de execução mostra chamadas incorretas da API OpenCL e problemas de acesso à memória.

Opencl

Para começar a usar o OpenCL, você precisa obter uma lista de plataformas disponíveis. As plataformas mais comuns são AMD, Intel e NVidia.


Em seguida, depois de escolher uma plataforma, você precisa selecionar o dispositivo de computação que essa plataforma representa:


E no final da fase preparatória, é necessário criar um contexto e filas dentro das quais a memória será alocada e os cálculos serão realizados. Por exemplo, um contexto que combina todos os dispositivos de computação de uma plataforma selecionada é criado da seguinte maneira:


Para baixar o código-fonte para um dispositivo de computação, ele deve ser compilado, a classe cl :: Program é destinada a isso.


Para descrever os parâmetros de uma função (kernel) executada em um dispositivo de computação, existe um modelo de cl: make_kernel.


Então tudo é simples: declaramos uma variável com o tipo de kernel, passamos o programa compilado e o nome do kernel computacional para ele, podemos iniciar o kernel quase como uma função normal.


O núcleo computacional do OpenCL é muito semelhante à opção 'openmp + block + optimization' da CPU, diferentemente da versão da CPU, o primeiro ciclo é controlado usando o OpenCL (o intervalo do ciclo é determinado pela variável global_range do código de inicialização do kernel, e o número da iteração atual está disponível no kernel usando a função get_global_id (0)). Primeiro, parte dos dados do corpo é carregada da memória local e processada. A memória local é comum a todos os encadeamentos do grupo, portanto, o download ocorre uma vez para o grupo e é processado por cada encadeamento no grupo e, desde a memória local é muito mais rápida que a global, os cálculos são muito mais rápidos.

Cuda

A implementação da plataforma NVidia CUDA é um pouco mais simples que o OpenCL, não precisamos criar o contexto do dispositivo e gerenciar a fila de execução (pelo menos até que desejemos fazer uma implementação com várias GPUs). Como no caso do OpenCL, precisamos alocar memória na GPU, copiar nossos dados nela e então podemos iniciar o núcleo da computação.

Você pode ler mais sobre como trabalhar com a CUDA aqui .


Diferentemente do OpenCL, o CUDA não especifica toda a gama de iterações (na implementação do OpenCL é global_range), mas define o tamanho da grade e os tamanhos dos blocos na grade; portanto, o cálculo do número do corpo atual muda ligeiramente; caso contrário, o kernel é muito semelhante ao OpenCL, com exceção de outros nomes. funções de sincronização e especificador para memória compartilhada. Outro recurso distintivo útil do CUDA é que podemos especificar o tamanho necessário da memória compartilhada quando o kernel é iniciado. Como na implementação do OpenCL, no início de cada bloco de iteração, copiamos parte dos dados na memória compartilhada e, em seguida, trabalhamos com essa memória de todas as threads do bloco.


Tabela 3. A dependência da função de tempo de cálculo (em segundos) $$$f_n$ no número de corpos interagindo $$$N$ para várias implementações de GPU e melhores implementações de CPU

$$$N$	4096	8192	16384	32768	65536	131072
openmp + block + otimização	0,0128	0,0495	0,194	0,774	---	---
OpenCL + meio NVidia K80	0,004	0,008	0,026	0,134	0,332	1,18
CUDA + meia NVidia K80	0,004	0,008	0,0245	0,111	0,291	1,13

Em geral, o resultado é decepcionante, apenas 5-6 vezes mais rápido que a implementação da CPU. Mesmo se realizarmos cálculos em todo o K80, obteremos aceleração até 12 vezes, mas desde Como a complexidade da tarefa é quadrática, em um tempo razoável, poderemos processar não 32768 corpos em interação, mas 131072, o que é apenas 4 vezes mais.

Função Aproximação $$$f_n$

Se você observar atentamente a função, definida pela força atrativa de dois corpos, fica claro que ela diminui quadraticamente com a distância. Portanto, podemos calcular com precisão a força de interação entre corpos próximos e aproximadamente entre corpos distantes. Uma abordagem famosa
é um algoritmo de código de árvore proposto por D. Barnes e P. Hat. Uma árvore octree é construída no espaço simulado, contendo em suas folhas as coordenadas e massas dos corpos simulados. Os nós pais contêm o centro de massa, a massa total dos nós filhos e o raio da esfera descrita ao redor dos corpos dos nós filhos. A raiz da árvore contém o centro de massa de todos os corpos, sua massa total e o raio da esfera descrita ao seu redor. Ao calcular a força de interação, a distância da raiz da árvore é considerada primeiro se a razão da distância do nó ao seu raio for maior do que alguma constante $$$\ lambda_ {crit}$ , então a raiz é considerada um corpo com coordenadas iguais às coordenadas do centro de massa dos corpos incluídos nela e uma massa igual à soma das massas dos nós filhos, se a razão for menor ou igual a $$$\ lambda_ {crit}$ , o procedimento é repetido recursivamente para cada nó filho. Se uma folha de uma árvore é alcançada, a força de interação é considerada da maneira usual. Assim, se um corpo está longe do grupo compacto de outros corpos, esse grupo é apresentado a ele como um corpo, e a força de interação é calculada não até cada corpo, mas apenas até um corpo. Devido a isso, a complexidade do algoritmo diminui com o $$$O (N ^ 2)$ antes $$$O (N \ log (N))$ ao custo de alguma perda de precisão.

Na minha abordagem, decidi não usar a árvore octree , mas a árvore kd, porque é mais fácil de usar e possui menor sobrecarga de armazenamento em comparação com o octree.

Vamos voltar à implementação na CPU. Um nó da árvore kd pode ser representado na forma de uma classe que contém ponteiros para os descendentes esquerdo e direito e informações sobre as coordenadas e a massa:


Com esse método de armazenamento da árvore, temos duas opções possíveis para percorrer a árvore: use recursão explícita ou use a pilha. Eu decidi pela segunda opção.


Como no caso da implementação "exata" da CPU, a função de cálculo de força é chamada para cada corpo. O ciclo em todos os corpos pode ser facilmente paralelizado usando diretivas OpenMP.

Mas as iterações de loop vizinhas, neste caso, se referem a partes completamente diferentes da árvore, o que não permite o uso eficiente do cache do processador. Para superar esse problema, todos os corpos precisam ser atravessados ​​não na ordem inicial, mas na ordem em que os corpos estão localizados nas folhas da árvore kd, iterações vizinhas ocorrerão para corpos que estão próximos no espaço e atravessarão a árvore por caminhos quase idênticos.


Essa implementação tem outro problema - não há uma maneira universal de paralelizar essa travessia de árvore. Mas podemos mudar completamente a maneira como a árvore é armazenada na memória - podemos armazenar todos os nós em uma matriz linear e abandonar completamente o armazenamento de ponteiros para os descendentes, por analogia com a construção de uma pilha binária . Ao iniciar a numeração de nós com $$$1$ se o nó estiver no número da célula $$$i$ , então o filho esquerdo dele está na cela $$$2i$ criança certa na célula $$$2i + 1$ e o pai na célula $$$i / 2$ . O nó direito correspondente à esquerda com o número $$$i$ terá um número $$$i + 1$ . A matriz em si terá um comprimento $$$2N$ e todos os nós das folhas estarão localizados no último $$$N$ Além disso, as células com espaçamento estreito serão localizadas nas células próximas da matriz. A função de passagem em árvore mudará um pouco, mas a estrutura permanece a mesma:


Mas essas não são todas as possibilidades que o armazenamento dos nós na matriz se abre para nós - podemos recusar a pilha ao rastrear. Para fazer isso, no ramo de código no qual vamos para os filhos do nó, adicionamos a função para calcular o próximo nó ( $$$next_ {up}$ ) e no ramo em que calculamos a força de interação, adicionamos o cálculo do próximo nó com a subárvore atual ignorada ( $$$skip_ {down}$ )

Para pular a subárvore atual, precisamos ir até a raiz (direção $$$D$ ), enquanto estamos no descendente direito, assim que alcançamos a esquerda, vamos para a subárvore direita correspondente (direção $$$R$ ), se chegarmos à raiz, a travessia da árvore será concluída.



Figura 8. Ignorando uma subárvore $$$skip_ {down}$ .

Para ir para o próximo nó, se possível, vá para o filho esquerdo (direção $$$U$ ) e, se não houver descendente, vá para o próximo nó 'de baixo' usando a função $$$skip_ {down}$ .



Figura 9. Transições para o próximo nó $$$next_ {up}$ .

Pode parecer que substituímos a recursão por um loop $$$enquanto$ em função $$$skip_ {down}$ , e essa substituição não fornece nada, mas vamos ver como determinar se um nó com um número $$$i$ descendente direito. Isso pode ser feito simplesmente verificando seu número ímpar (o filho certo está na célula $$$2i + 1$ ), para isso basta calcular $$$i \ e 1$ . I.e. nós dividimos o número $$$i$ dois se o bit menos significativo estiver definido como um. Mas isso pode ser feito sem um loop, em muitos processadores existe uma instrução find first set que retorna a posição do primeiro bit configurado, usando-o, transformamos o loop em quatro instruções do processador.


Depois disso, podemos excluir a pilha da função de deslocamento em árvore e substituí-la por um par $$$skip_ {down} + próximo_ {up}$ , depois disso a função parece ainda mais simples:


No total, obtivemos seis combinações de deslocamento de árvore e cálculo de força. Compare essas abordagens em termos de tempo e qualidade de computação. Tomamos como medida de qualidade a mudança relativa na energia total do sistema após 100 iterações. Como sistema modelo, tomamos duas "galáxias" interagindo consistindo em$$$16384$corpos cada.

Tabela 4. Combinações do método transversal de árvore e cálculo de força

Tipo de cálculo de deslocamento / força em árvore	Árvore com pilha	'Pilha' com uma pilha	'Heap' sem uma pilha
Iterações por número do corpo	ciclo + árvore	ciclo + pilha	cycle + heapstackless
Bypass de folhas	nestedtree + tree	nestedtree + heap	nestedtree + heapstackless

a) Dependência do tempo de cálculo da função $$$f_n$ da razão entre a distância crítica e o nó da árvore e seu raio ( $$$\lambda_{crit}$ )	b) A dependência da mudança relativa de energia no sistema na razão entre a distância crítica e o nó da árvore e seu raio ( $$$\lambda_{crit}$ )
Figura 10. Resultados do cálculo da função $$$f_n$ de várias maneiras na CPU (número de corpos $$$N=32768$ )

Pode-se ver que a implementação de 'nestedtree + tree' está irremediavelmente atrasada, porque falta concorrência. Implementações com a localização de nós de árvore na matriz e a indexação como em um heap binário estão na liderança. A mudança relativa de energia é insignificante para todas as variantes com$$$\lambda_{crit} > 1$ .Também na fig. 10a mostra que para ($$$\lambda_{crit} < 10$ ) todas as variantes do cálculo da função $$$f_n$ ultrapassar significativamente a velocidade da versão mais otimizada do cálculo exato ('openmp + bloco + otimização'), com um aumento adicional $$$\lambda_{crit}$ implementações com uma árvore perdem a versão exata.

Travessal em árvore da GPU

Tentei ignorar a árvore na GPU usando a tecnologia OpenCL e CUDA. A opção de armazenar nós na forma de uma árvore foi imediatamente descartada e apenas as opções foram deixadas com o armazenamento da árvore em uma matriz com indexação como em uma pilha binária. Em geral, a implementação do núcleo de computação não é muito diferente da versão da CPU.


Na primeira versão, o percurso da árvore começou na ordem de numeração dos corpos na matriz original; portanto, os threads vizinhos ignoraram partes completamente diferentes da árvore, o que afetou negativamente o desempenho do cache de memória da GPU. Portanto, na segunda modalidade, uma travessia foi aplicada, começando do topo da árvore, neste caso, os segmentos vizinhos começam a atravessar a árvore pelo mesmo caminho, porque as copas das árvores vizinhas estão próximas e no espaço. Também é importante escolhermos a numeração na matriz de nós da árvore não do zero, mas de uma, neste caso, as folhas da árvore são armazenadas na segunda metade da matriz e, com o número de corpos igual à potência de dois, teremos acesso igual à memória pelo índice tn1.


Ao percorrer a ordem de numeração dos nós da árvore, obtivemos um aumento de desempenho. Mas esta opção também pode ser melhorada. A memória global na qual os nós da árvore estão localizados atualmente é otimizada para acesso compartilhado , ou seja, threads de um grupo devem ler palavras localizadas em um bloco de memória. No nosso caso, o percurso da árvore começa nos mesmos caminhos, e solicitamos os mesmos dados com todos os threads do grupo, mas quanto mais nos aprofundamos na árvore, os caminhos dos threads vizinhos divergem cada vez mais, e precisamos solicitar dados diferentes, o que reduz o desempenho do subsistema várias vezes. memória. Mas os nós de cada subárvore pertencente ao mesmo nível estão localizados em células de memória relativamente próximas. I.e.ao percorrer o restante da árvore, os segmentos vizinhos do núcleo de computação não acessam os mesmos nós da árvore, mas localizados na memória. Para otimizar esse acesso à memória, pode ser usada memória de textura. Mas há um problema. No momento, as texturas não suportam o trabalho com dados de precisão dupla (queremos calcular com precisão). Mas no CUDA há uma função __hiloint2double que coleta um número de precisão dupla de dois números inteiros.


Nesse caso, duas implementações foram feitas, em uma cada elemento da árvore (x, y, z, tree_crit_r2) foi solicitado independentemente, e na segunda implementação essas solicitações foram combinadas. A solicitação para a massa do nó ocorre com muito menos frequência, apenas se a condição r2> tree_crit_r2 [curr] for atendida , portanto, não faz sentido combinar essa solicitação com as outras. Outro recurso útil da estrutura CUDA é a capacidade de controlar a proporção dos tamanhos do cache L1 e o tamanho da memória compartilhada ( cudaFuncSetCacheConfig ). No caso de passagem de árvore, não usamos memória compartilhada, portanto, podemos aumentar o cache L1 em ​​detrimento dele.


Uma análise do programa no nvprof profiler mostrou que, mesmo ao usar a memória de textura para armazenar uma árvore, ainda há uma carga muito alta na memória global.

De fato, no CUDA, toda a memória do kernel endereçada a endereços 'calculados' é armazenada na memória global e, portanto, a pilha necessária para percorrer a árvore está localizada na memória global e consome uma parte significativa da largura de banda do chip de memória, porque a pilha possui cada thread de execução e há muitos threads.

Mas, felizmente, já sabemos como atravessar uma árvore sem usar uma pilha. Complementando o núcleo computacional anterior com as funções de calcular o próximo nó da árvore, obtemos um novo kernel, além disso, mais compacto.


O desempenho dos núcleos em execução na GPU depende muito do tamanho dos blocos nos quais dividimos a tarefa. Depende desse tamanho, quantos registros, memória local e outros recursos estarão disponíveis para cada thread de computação. Você também deve ter em mente que, enquanto aguarda o acesso à memória em um thread, outro thread pode executar cálculos no processador de sombreador, portanto, com um número suficiente de threads em execução simultaneamente em um processador, o tempo de acesso à memória ficará oculto por trás dos cálculos. Portanto, antes de comparar o desempenho de nossos núcleos, precisamos calcular o tamanho ideal do bloco para cada um deles. Vamos fazer uma comparação na metade disponível da NVidia K80 para nós.

Tabela 5. A dependência da função de tempo de cálculo (em segundos) $$$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

Uma situação difícil, mas, diferentemente da versão da CPU do percurso da árvore, é claro que cada etapa de otimização traz resultados tangíveis. A implementação do 'opencl + heap + cycle' é quase 6 vezes mais lenta que uma solução exata com cálculo de função completa$$$f_n$ .Uma implementação de 'opencl + heap + nested', na qual uma travessia de árvore inicia nos nós vizinhos, 1,4 vezes mais rápido que o anterior, porque A memória cache é usada melhor. Na implementação de 'cuda + heap + anested', o cache L1 foi aumentado em detrimento da memória compartilhada, o que aumentou a velocidade em 1,4 vezes, embora seja possível que na implementação do cuda o kernel de compilação seja mais otimizado (nas versões de 'opencl + denso' e 'cuda os núcleos + densos são idênticos e o desempenho da versão do cuda é ~ 1,2 vezes maior). O aumento mais significativo na velocidade computacional (3,8 vezes) é alcançado quando a árvore está localizada na memória de textura e as consultas aos elementos do nó da árvore são combinadas. A implementação com o deslocamento em árvore sem usar a pilha 'cuda + heap + nested + tex + stackless' também é 1,4 vezes mais rápida.nele, toda a largura de banda do barramento de memória é usada apenas para acessar dados sobre nós da árvore e não é gasta na pilha. Assim, com$$$\lambda_{crit}=10$ Foi possível alcançar a aceleração duas vezes em comparação com o cálculo completo da função $$$f_n$ . Mas $$$\lambda_{crit}=10$ um valor excessivamente grande do parâmetro, no gráfico da mudança relativa de energia no sistema, da razão entre a distância crítica e o nó da árvore e seu raio para a implementação da CPU, pode ser visto que valores mais baixos podem ser usados $$$\lambda_{crit}$sem uma perda visível de precisão. Vamos tentar variar$$$\lambda_{crit}$ com os tamanhos de bloco ideais que determinamos na etapa anterior.

a) $$$N=128K$	b) $$$N=1M$
Fig. 11. Dependência do tempo de cálculo da função. $$$f_n$ da razão entre a distância crítica e o nó da árvore e seu raio ( $$$\lambda_{crit}$ ) para várias implementações de passeio em árvore da GPU

Pode-se ver que, para pequenas $$$\lambda_{crit}$ todos os métodos de cálculo de funções $$$f_n$vá para valores próximos, determinados no momento da construção da árvore kd e da preparação dos dados para a GPU. Além disso, o tempo de construção de uma árvore contribui significativamente para o tempo total de$$$\lambda_{crit}\leq 4$, esse tempo pode ser negligenciado. É interessante notar que quando$$$N=128K$ o desempenho melhora novamente ao atingir $$$\lambda_{crit}=1024$provavelmente isso se deve ao fato de que todos os threads da GPU percorrem os mesmos vértices da árvore ao mesmo tempo, o que melhora o uso do cache L1 e elimina completamente a 'divergência de ramificação' . Também é visto que, com o melhor desempenho para uma implementação sem uma pilha (cuda + heap + nested + tex + stackless), ele supera a versão com a pilha em aproximadamente$$$1.4-1.5$vezes. Outras implementações são várias vezes mais lentas. Para consolidar o resultado, calcularemos o tempo em uma GPU com uma arquitetura mais nova.


Ao usar o GeForce GTX 1080 Ti para cálculo, a diferença entre as caminhadas pela árvore com uma pilha e sem uma pilha é de até duas vezes, desde que negligenciemos o tempo necessário para construir a árvore. Esse fato nos leva a garantir que, no futuro, portemos para a GPU e construamos uma árvore. Uma comparação das tabelas 5-7 mostra que não há um tamanho ideal de bloco ideal para um número diferente de corpos e, ainda mais, para diferentes arquiteturas de GPU, além disso, a diferença de tempo de cálculo pode atingir várias vezes, mesmo se você não considerar os valores de limite. Portanto, antes de cálculos demorados, faz sentido determinar o tamanho ideal do bloco para cada configuração.

O principal que alcançamos é a capacidade de modelar a interação pareada de pouco mais de um milhão de corpos ($$$2^{20}$) por um tempo razoável em um, não na GPU mais recente. Em GPUs mais recentes (Tesla V100), aparentemente, será possível processar já cerca de dois milhões de corpos em interação em um tempo razoável, porque seu desempenho é cerca de quatro vezes maior que a metade do Tesla K80. Embora esse número também seja incomparável com o número de estrelas, mesmo em galáxias anãs, como a Pequena Nuvem de Magalhães , mas, na minha opinião, é impressionante.

Conclusão

O uso de métodos incorporados para resolver equações diferenciais torna possível resolver problemas semelhantes com um alto nível de precisão a um custo de tempo relativamente pequeno, e a aplicação de algoritmos para aproximar a função de calcular as forças de atração aos pares nos permite processar volumes colossais de corpos em interação. Assim, diferentemente dos alienígenas das "Tarefas de três corpos", somos capazes de resolver o problema$$$N$ corpos e para um pequeno número de corpos e aglomerados de estrelas inteiros, embora ao custo de alguma perda de precisão.

Visualização

Para quem leu até o fim, darei mais alguns vídeos com visualização da evolução dos sistemas dos corpos.

Simulação de uma colisão de duas galáxias. O número total de corpos é de 60 mil.


Modelando a evolução de uma galáxia que consiste em um milhão de estrelas. No centro está um corpo com uma massa de 99% do total. Partículas únicas são quase indistinguíveis. Já é mais como uma onda em uma gota de líquido. Colorido de acordo com o módulo de velocidade. Baixa velocidade - azul, verde - médio, vermelho - alto. Vê-se que no centro a velocidade é mais alta e diminui gradualmente para as bordas, e a mais baixa no plano equatorial.


Uma simulação mais "dinâmica" da evolução de uma galáxia de milhões de estrelas. No centro está um corpo com uma massa de 10% do total. O corpo central afeta o resto significativamente menos. Primeiro, as "estrelas" se separam, depois se reúnem em vários aglomerados e, no final, formam novamente um grande aglomerado.

No processo de modelagem, cerca de 5% das "estrelas" deixaram a região inicial irrevogavelmente.


No décimo segundo, assemelha-se muito à galáxia de roda dentada existente .

O código do projeto pode ser encontrado no github .