DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE ESTRUTURAS DE VOXELS USANDO UNIDADES DE PROCESSAMENTO GRÁFICO - PDF
DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE ESTRUTURAS DE VOXELS USANDO UNIDADES DE PROCESSAMENTO GRÁFICO
Download "DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE ESTRUTURAS DE VOXELS USANDO UNIDADES DE PROCESSAMENTO GRÁFICO"
Lucca Desconhecida Meneses
1 AUTARQUIA ASSOCIADA À UNIVERSIDADE DE SÃO PAULO DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE FÓTONS EM ESTRUTURAS DE VOXELS USANDO UNIDADES DE PROCESSAMENTO GRÁFICO Murillo Bellezzo Dissertação apresentada como parte dos requisitos para obtenção do Grau de Mestre em Ciências na Área de Tecnologia Nuclear - Reatores Orientador: Prof. Dr. Hélio Yoriyaz São Paulo 2014
2 INSTITUTO DE PESQUISAS ENERGÉTICAS E NUCLEARES AUTARQUIA ASSOCIADA À UNIVERSIDADE DE SÃO PAULO DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE FÓTONS EM ESTRUTURAS DE VOXELS USANDO UNIDADES DE PROCESSAMENTO GRÁFICO Murillo Bellezzo Dissertação apresentada como parte dos requisitos para obtenção do Grau de Mestre em Ciências na Área de Tecnologia Nuclear - Reatores Orientador: Prof. Dr. Hélio Yoriyaz São Paulo 2014
3 Resumo Sendo o método mais preciso para estimar a dose absorvida em radioterapia, o Método de Monte Carlo (MMC) tem sido amplamente utilizado no planejamento de tratamento radioterápico. No entanto, a sua eficiência pode ser melhorada para aplicações clínicas de rotina. Nesta dissertação é apresentado o código CUBMC, um código de Monte Carlo que simula o transporte de fótons para cálculo de dose, desenvolvido na plataforma CUDA (Compute Unified Device Architecture). A simulação de eventos físicos é baseada no algoritmo presente no código PENELOPE, e as tabelas de seção de choque utilizadas são geradas pela rotina MATERIAL, também presente no código PENELOPE. Os fótons são transportados em objetos simuladores descritos por voxels. Existem duas abordagens distintas utilizadas para a simulação. A primeira delas obriga o fóton a realizar uma parada toda vez que cruza a fronteira de um voxel, a segunda é pelo Método de Woodcock, onde o fóton ignora a existência de fronteiras e é transportado em um meio homogêneo fictício. O código CUBMC tem como objetivo ser uma opção de código simulador que, ao utilizar a capacidade de processamento paralelo de unidades de processamento gráfico (GPU), apresente alto desempenho em máquinas compactas e de baixo custo, podendo assim ser aplicado em casos clínicos e incorporado a sistemas de planejamento de tratamento em radioterapia. Palavras-chaves: GPU; dosimetria; simulação de Monte Carlo; radioterapia; transporte de fótons; PENELOPE; MCNP5
4 Abstract As the most accurate method to estimate absorbed dose in radiotherapy, Monte Carlo Method (MCM) has been widely used in radiotherapy treatment planning. Nevertheless, its efficiency can be improved for clinical routine applications. In this thesis, the CUBMC code is presented, a GPU-based MC photon transport algorithm for dose calculation under the Compute Unified Device Architecture (CUDA) platform. The simulation of physical events is based on the algorithm used in PENELOPE, and the cross section table used is the one generated by the MATERIAL routine, also present in PENE- LOPE code. Photons are transported in voxel-based geometries with different compositions. There are two distinct approaches used for transport simulation. The first of them forces the photon to stop at every voxel frontier, the second one is the Woodcock method, where the photon ignores the existence of borders and travels in homogeneous fictitious media. The CUBMC code aims to be an alternative of Monte Carlo simulator code that, by using the capability of parallel processing of graphics processing units (GPU), provide high performance simulations in low cost compact machines, and thus can be applied in clinical cases and incorporated in treatment planning systems for radiotherapy. Keywords: GPU; dosimetry; Monte Carlo simulation; radiotherapy; photon transport; PENELOPE; MCNP5
5 Sumário 1 Introdução 1 2 Objetivos 3 3 Revisão Bibliográfica 4 4 Materiais e Métodos Unidade de Processamento Gráfico (GPU) CUDA (Compute Unified Device Architecture) Arquitetura de uma GPU Tipos de Memória Ocupação Método de Monte Carlo Técnicas de amostragem Método da Inversão ou Método Direto Distribuições Discretas Técnica da rejeição Gerador de números pseudoaleatórios RANECU Paralelização do gerador de números pseudoaleatórios Transporte de fótons Sorteio do passo Propagação da partícula
6 Sorteio da interação Cálculo dos tallies AMIGOBrachy Resultados Gerador de Números Pseudoaleatórios (PRNG) Comparativo de desempenho entre o RANECU e sua versão paralelizada Função Rand Cálculo de π CUBMC Aspectos Gerais Geometria Tabelas de Seção de choque Modelos de fonte de radiação Sorteio do passo Fronteiras Interações Estimadores (Tallies) Finalização da história da partícula Método de Woodcock Simulações Discussões 71 7 Conclusão 73 Referências Bibliográficas 76
7 1 Introdução O Método de Monte Carlo (MMC) é um método estatístico baseado no uso de números aleatórios que pode ser aplicado para solucionar qualquer problema cuja representação possa ser descrita por uma sucessão de processos estocásticos. Ele é amplamente utilizado na simulação do transporte de radiação na matéria para a obtenção de grandezas físicas de interesse, tais como fluxo de partículas e energia depositada. Devido à natureza estatística de suas simulações, o Método de Monte Carlo pode ser facilmente paralelizado. Dentre os programas computacionais existentes baseados no Método de Monte Carlo destacam-se o MCNP6 [11], PENELOPE [5] e GEANT4 [9]. Estes programas foram desenvolvidos e aperfeiçoados ao longo de décadas de pesquisas na área e hoje são amplamente usados na área de física médica. Entretanto, são ferramentas de uso geral, muito complexas e de difícil adaptação a novas plataformas, exigindo uma demanda de estudos e experiência em simulação por parte dos usuários. Esforços vem sendo realizados para a adaptação destes códigos ao processamento paralelo como, por exemplo, no caso do programa PENELOPE [4]. Outros programas como o MCNP já estão adaptados para serem executados paralelamente, seja em computadores com processadores de multiplos núcleos ou em cluster s, que são supercomputadores compostos por centenas ou milhares de processadores que se comunicam de modo a realizar processamento paralelo. Cluster s apresentam um alto poder de processamento, no entanto, são pouco acessíveis devido ao seu alto custo. Uma alternativa aos clusters são as unidades de processamento gráfico (GPUs), que possuem uma arquitetura projetada para o processamento paralelo. Tendo sua evolução alavancada pela indústria de jogos, a capacidade de processamento aritmético das GPUs ultrapassou a das unidades de processamento central (CPUs) e softwares mais recentes fazem uso dessa capacidade de processamento aritmético em paralelo para ganho de desempenho em seus processos, como podemos observar na figura 1.1. Visando tirar proveito dessa alta capacidade de processamento aritmético paralelo das GPUs, surgiram novos códigos de programação capazes de direcionar parte do
8 1 Introdução 2 Figura 1.1: Comparativo entre operações de pontos flutuantes por segundo entre as principais CPUs e GPUs fluxo de suas operações para as GPUs. Dentre eles, se destaca o CUDA (Compute Unified Device Architecture), que é a linguagem usada neste trabalho. O presente trabalho consiste no desenvolvimento do código CUBMC (CUDA Based Monte Carlo), um software que usa o poder de execução em paralelo das unidades de processamento gráfico para acelerar simulações computacionais que utilizam o Método de Monte Carlo para o transporte de fótons.
9 2 Objetivos O objetivo deste trabalho é o desenvolvimento do software CUBMC (CUDA Based Monte Carlo), um código de Monte Carlo destinado ao transporte de fótons em objetos simuladores descritos por voxels, para que ele seja uma opção de código simulador que, ao utilizar a capacidade de processamento paralelo de unidades de processamento gráfico (GPUs), apresente alto desempenho em máquinas compactas e de baixo custo, podendo assim ser aplicado em casos clínicos e incorporado a sistemas de planejamento de tratamento de radioterapia. Dessa forma, a opção de trabalhar com objetos simuladores definidos por voxels busca facilitar a descrição de modelos anatômicos provenientes de imagens de Tomografia Computadorizada (TC) ou Ressonância Magnética (RM), que são facilmente extrapoladas para voxels. Visando obter ganho de desempenho, o CUBMC destina-se exclusivamente à simulação do transporte de fótons, as partículas secundárias geradas nas interações não são transportadas e sua energia é depositada localmente. A física envolvendo a simulação do transporte da radiação presente no código foi simplificada ao máximo dentro do limite em que ainda apresente um resultado adequado no intervalo energético de interesse para os tratamentos de radioterapia e braquiterapia (1keV 9MeV). Dentre as interações simuladas (Efeito Fotoelétrico, Espalhamento Rayleigh, Espalhamento Compton e Produção de Pares), apenas o Espalhamento Compton apresenta um algoritmo de simulação mais abrangente, ele foi extraído sem alterações do código PENELOPE. As outras interações são extremamente simplificadas e serão implementadas à medida que o código CUBMC for aprimorado. As tabelas de seção de choque são geradas pela rotina MATERIAL, também pertencente ao código PENELOPE. O código CUBMC permite que o transporte dos fótons seja simulado segundo dois métodos. O método convencional em que a partícula é obrigada a realizar uma parada toda vez que encontrar um fronteira da geometria (chamado neste trabalho de método de paradas ), e o método em que a partícula é transportada em um meio homogêneo fictício, ignorando as fronteiras da geometria (Método de Woodcock).
10 3 Revisão Bibliográfica As simulações realizadas pelo método de Monte Carlo tem sido amplamente usadas na área de física médica nos últimos 60 anos. Se pesquisarmos pelo termo Monte Carlo na PubMed (http://www.ncbi.nlm.nih.gov/pubmed) em fevereiro de 2014, teremos mais de 36 mil retornos, sendo que o primeiro artigo data de setembro de 1949 [21]. Entre os anos de 1967 e 2000, o número de publicações contendo o termo Monte Carlo em seu título ou resumo presentes nas revistas Physics in Medicine and Biology (PMB) e Medical Physics dobraram a cada, aproximadamente, 5 anos [24], conforme ilustra a figura 3.1. Figura 3.1: Número de publicações contendo o termo Monte Carlo em seu título ou resumo presentes nas revistas Physics in Medicine and Biology (PMB) e Medical Physics ao longo dos anos [24] Esse crescimento pode ser atribuído ao aumento da capacidade de processamento e redução de custo dos computadores, que alavancou o desenvolvimento de softwares envolvendo o Método de Monte Carlo para simular o transporte de radiação. O emprego destes softwares pode ser destacado para o transporte de radiação em cálculos dosimétricos
11 3 Revisão Bibliográfica 5 de braquiterapia e radioterapia, sendo parte integrante fundamental em alguns sistemas de planejamento de tratamento [2]. Para que o tempo de simulação dos códigos de Monte Carlo fossem viáveis para aplicação clínica em sistemas de planejamento de tratamento, foram desenvolvidos métodos para que esses códigos fossem executados mais rapidamente. Dentre eles, podemos destacar o VMC++ (Voxel Monte Carlo) [15] que é semelhante ao EGSnrc [16], e resultou em simulações ao mínimo uma ordem de grandeza mais rápidas [24]. Também temos o DPM [26], que teve os modelos físicos para as simulações simplificados para que a execução fosse mais rápida, mas ainda obtivesse resultados corretos quando aplicado a sistemas de planejamento radioterápicos e o MCDOSE [20]. Ainda assim, nenhum deles obteve a mesma eficiência que o código VMC++. A descrição da geometria por voxels mostrou-se adequada para objetos simuladores provenientes de modelos anatômicos, uma vez que imagens de tomografia computadorizada (CT) e ressonância magnética (RM) são facilmente exportadas para modelos tridimensionais voxelizados, além de permitir descrições mais detalhadas do corpo humano [18] [17] e específicas para cada paciente [8]. Outra forma de se acelerar uma simulação de Monte Carlo é paralelizar o processo. Uma vez que as partículas simuladas são independentes entre si, elas podem ser simuladas simultaneamente por núcleos processadores distintos. Esforços foram feitos para paralelizar códigos de Monte Carlo já existentes, como foi feito com o PENELOPE [10], ou o desenvolvimento de novas versões dos códigos englobando a possibilidade de processamento paralelo como uma de suas atualizações, como o MCNP5 [7]. A execução em paralelo possibilitou que códigos de Monte Carlo mais abrangentes fossem executados em tempos comparáveis aos de códigos simplificados. No entanto, existia a limitação financeira devido ao alto custo de computadores com grande capacidade de processamento paralelo. Com o advento do uso das Unidades de Processamento Gráfico (GPUs) 1 como processadores auxiliares, o desenvolvimento ou adaptação de softwares de Monte Carlo otimizados para processamento paralelo foi alavancado [23]. 1 GPUs apresentam uma elevada capacidade de processamento paralelo a um baixo custo
12 3 Revisão Bibliográfica 6 Códigos como o PENELOPE [3] e DPM [13] foram migrados para a execução em GPU, contudo, os ganhos de desempenho foram modestos quando comparados ao potencial de ganho. Isso devido, principalmente, a serialização que ocorre ao se gerar partículas secundárias nesses códigos. Já o código GPUMCD [12] obteve ganho de desempenho de 200x frente ao DPM após corrigir o tratamento das partículas secundárias. Para simulações que consideram apenas o transporte dos fótons, ganhos de desempenho de até 1000x foram obtidos para geometrias simples [1]. Em 2011, o código gctb [14], desenvolvido sob a plataforma CUDA [22], obteve ganhos de 76.6x frente ao código EGSnrc para simulações para o modelo anatômico de cabeça Zubal [28], e de 400x para uma geometria composta por um meio homogêneo de água. Embora o desenvolvimento de códigos de Monte Carlo otimizados para execução em GPU ainda seja incipiente, ganhos expressivos já foram obtidos. O poder de processamento paralelo das GPUs ainda é crescente, e seu baixo custo faz com que exista um grande potencial de que elas venham a substituir os cluster s para simulações de Monte Carlo de alta performance.
13 4 Materiais e Métodos 4.1 Unidade de Processamento Gráfico (GPU) A GPU surgiu como um processador auxiliar destinado para processamento gráfico, com o intuito de reduzir a quantidade de processamento lógico aritmético realizado pela unidade de processamento central (CPU). Analisando a arquitetura de uma CPU e de uma GPU (fig. 4.1), podemos observar que uma CPU tem a sua área de controle desenvolvida enquanto suas unidades lógicas aritméticas (ALU) são escassas. Já a GPU tem apenas alguns poucos núcleos de controle e muitas unidades ALU. Figura 4.1: Diferença arquitetônica entre CPU e GPU Tal diferença na arquitetura faz com que os processadores centrais sejam dispositivos otimizados para o gerenciamento paralelo de tarefas mas pouco eficientes para realizar operações aritméticas. Já as unidades de processamento gráfico possuem pouca capacidade no que tange ao processamento de múltiplas tarefas, mas são altamente eficientes na execução de operações aritméticas. Visando tirar proveito dessa alta capacidade de processamento aritmético paralelo das GPUs, surgiram novos códigos de programação capazes de direcionar parte do fluxo de suas operações para as GPUs. Dentre eles, se destaca o CUDA (Compute Unified Device Architecture), que é a linguagem usada nesse trabalho.
14 4.1 Unidade de Processamento Gráfico (GPU) CUDA (Compute Unified Device Architecture) CUDA é uma linguagem de programação desenvolvida pela NVIDIA voltada para as GPUs da própria empresa. Visando uma rápida curva de aprendizagem para programadores já familiarizados com outros códigos, CUDA é uma plataforma com versões que compreendem os comandos e as bibliotecas existentes em C, C++, C# e Fortran. [22] CUDA C é uma linguagem semelhante à linguagem C/C++, com o diferencial de possuir ferramentas com as quais é possível estabelecer uma comunicação entre a GPU e a CPU. Dessa forma, as partes do código que são responsáveis por operações aritméticas podem ser executadas na GPU. Entretanto, para um desempenho otimizado, o código deve ter seus parâmetros (discutidos em detalhes mais adiante) ajustados segundo a capacidade da GPU na qual será executado. Diferentes gerações de GPUs NVIDIA apresentam diferentes compatibilidades CUDA. Essas distinções fazem com que a forma como a GPU é capaz de gerenciar os processos em paralelo, e a memória disponível para cada processo seja diferente Arquitetura de uma GPU Embora possuam a mesma estrutura geral, a arquitetura das unidades de processamento gráfico compatíveis com a tecnologia CUDA varia conforme a geração das GPUs. Cada GPU compatível com a tecnologia CUDA é dividida fisicamente por um grupo de multiprocessadores compostos por múltiplos núcleos, de modo que, uma GPU tenha de centenas à milhares de núcleos CUDA. Essa divisão é intrínseca de cada GPU e não pode ser alterada pelo código. Quando uma GPU vai executar alguma tarefa, a execução se dá através de blocos chamados warps compostos por 32 thread s (menor unidade de processamento de uma GPU). As GPUs já são projetadas para processamento paralelo, e mesmo que ela tenha que processar apenas um grupo de informações, ela vai realizar o mesmo procedimento que seria necessário para processar 32 grupos. Para se ter um exemplo na figura 4.2 temos o esquema da GPU GTX 560 Ti
15 4.1 Unidade de Processamento Gráfico (GPU) 9 (usada nesse trabalho). Ela possui 8 multiprocessadores de 48 núcleos cada, totalizando 384 núcleos. Cada núcleo pode executar um warp (de 32 threads) por vez. Assim, a ocupação máxima da GPU é de threads. Figura 4.2: Modelo esquemático da arquitetura de uma GPU GTX 560 Ti (à esquerda) e de um de seus multiprocessadores (à direita) Além da divisão física presente nas GPUs, existem parâmetros virtuais que definem uma hierarquia referente à execução de alguma função (kernel) enviada pela CPU (host). A menor unidade de processamento ao se executar um kernel é o thread. Osthread s não são núcleos físicos, são organizados em blocos, e são responsáveis por realizar as operações aritméticas de uma determinada função. Sua organização pode ser unidimensional, bidimensional e até tridimensional. Diferente formas de agrupamento podem facilitar a identificação dos thread s para um determinado problema. Agrupamentos bidimensionais, por exemplo, podem ser interessantes ao se trabalhar com matrizes, no entanto, a organização unidimensional tende a promover um acesso à memória de maneira mais eficiente. Todos os threads de um mesmo bloco tem acesso a uma memória comum (memória compartilhada ou shared ) e devem ser executadas pelo mesmo kernel.
16 4.1 Unidade de Processamento Gráfico (GPU) 10 Os blocos são organizados em grades (figura 4.3) sendo que, os blocos de uma mesma grade também são executados pelo mesmo kernel, mas são independentes entre si no que diz respeito à memória e não são executados em sincronia. Desse modo, os threads de um bloco não cooperam com threads de outro bloco. Cada núcleo CUDA será responsável pela execução de um bloco. Sendo assim, existe um limite para o número de blocos que um multiprocessador pode executar simultaneamente. Caso o número de blocos a serem executados exceda o número de núcleos existentes na GPU, é criada uma fila e os blocos serão executados à medida que outros blocos encerrem suas atividades. Figura 4.3: Esquema das divisões hierárquicas na execução de um kernel Tipos de Memória Para um melhor aproveitamento da massiva capacidade de paralelização das unidades de processamento gráfico é necessário usar de forma adequada os diferentes tipos de memória nelas presentes, pois cada tipo de memória apresenta características diferentes,
17 4.1 Unidade de Processamento Gráfico (GPU) 11 como forma de escrita e leitura, capacidade de armazenamento, tempo necessário para a leitura e escrita e otimização para acesso simultâneo. A figura 4.4 ilustra onde cada tipo de memória fica alocada fisicamente na GPU e quem pode acessá-la. Figura 4.4: Modelo esquemático das diferentes memórias presentes na GPU Memória Global A memória global é a única memória que pode ser acessada para leitura e escrita, tanto pela CPU (host) quanto pela GPU (device). Desse modo, torna-se a principal responsável pela troca de dados entre esses processadores. Apresenta a maior capacidade de alocação da GPU (Gb), e a informação nela contida persiste durante a execução inteira do programa, podendo ser acessada por todos os elementos da grade. No entanto, não é otimizada para acessos simultâneos e exige um tempo longo para ser lida ( ciclos).
18 4.1 Unidade de Processamento Gráfico (GPU) 12 A memória global é a memória usada como padrão para a comunicação entre CPU e GPU caso o código não especifique o uso de outra memória. Memória Constante Semelhante à memória global, pode ser lida e escrita pelo host mas a GPU só tem permissão para leitura. Assim como na memória global, a informação contida na memória constante persiste durante a execução inteira do aplicativo. Dessa forma, é utilizada para o armazenamento de informações que não serão alteradas durante a execução do kernel, podendo ser acessada por todos os elementos da grade. É otimizada para acessos simultâneos e possui um tempo de leitura moderado ( ciclos), contudo, apresenta capacidade limitada de armazenamento (64Kb/multiprocessador). Memória de Registro Memória usada automaticamente pelo thread para alocação de variáveis locais. Possui o menor tempo de acesso ( 1 ciclo) mas é visível apenas pelo thread. Cada multiprocessador tem uma quantidade fixa de memória de registro, que é dividida entre os threads nele presentes. Tem a duração da execução do kernel. O uso correto da memória de registro pode acarretar em uma significativa redução no tempo de execução do código. Memória Compartilhada (Shared) Semelhante à memória de registro, mas pode ser acessada por todos os threads de um mesmo bloco. Embora tenha acesso mais lento que o da memória de registro (de 3 a 6 vezes mais lento), é especialmente útil para a comunicação entre os threads sem que seja necessário acessar a memória global (muito mais lenta). Memória de Textura A memória de textura é a própria memória global, no entanto, é acessada por um cache dedicado apenas à leitura e otimizado para um acesso mais rápido de elementos vizinhos (no sistema de coordenadas definido na memória de textura). A memória de textura é uma ótima opção para se armazenar as informações sobre o material, por exemplo, uma vez que voxels vizinhos tendem a ser acessados com
19 4.2 Método de Monte Carlo 13 maior frequência durante a história de uma partícula Ocupação Ocupação é definida como a razão entre o número de warps sendo processados e o número máximo de warps permitido pela GPU. A princípio, uma maior ocupação significaria um maior uso dos recursos da GPU. Entretanto, isso não implica necessariamente em um melhor desempenho. Cada multiprocessador possui um limite para memória de registro e para memória compartilhada (que são as memórias de acesso mais rápido da GPU). Sendo assim, quanto menos thread s forem executados simultaneamente por processador, maior a quantidade de memória disponível para cada um deles. Sacrificar o grau de paralelização em detrimento de mais registros disponíveis por thread pode ser lucrativo em relação à eficiência do código devido ao acesso à memória de registro ser muito mais rápido que o acesso à memória global. Contudo, também existe um limite máximo de memória por thread. Tais peculiaridades fazem com que existam planilhas para o cálculo da ocupação da GPU dada uma configuração. No entanto, testes empíricos acabam sendo necessários para se obter a configuração de melhor desempenho para um determinado problema. 4.2 Método de Monte Carlo O Método de Monte Carlo (MMC) trata-se de uma classe de métodos numéricos baseados no uso de números aleatórios. Recebeu esse nome na década de 1940 por cientistas integrantes do Projeto Manhattan, em Los Alamos, devido à aleatoriedade em suas simulações estatísticas, também presente nos jogos de azar, e por causa do bairro de Monte Carlo, pertencente ao principado de Mônaco, ser conhecido na época como a capital mundial dos jogos de azar. Desde então, o método vem sendo amplamente usado para solucionar problemas complexos de física e matemática sendo, muitas vezes, o único meio viável para problemas de múltiplas variáveis independentes para os quais os métodos numéricos con-
20 4.2 Método de Monte Carlo 14 vencionais exigiriam um tempo computacional, e/ou quantidade de memória demasiadamente elevados. O Método de Monte Carlo é um método estatístico que pode ser aplicado para solucionar qualquer problema cuja representação possa ser descrita por uma sucessão de processos estocásticos. Em muitas aplicações o processo físico pode ser simulado sem a necessidade de se solucionar as equações matemáticas que o descrevem, sendo apenas necessário que ele seja representado por funções de densidade de probabilidade (PDFs). A amostragem aleatória dessas PDFs simula de forma determinística uma sequência específica de eventos e o resultado dessa simulação passa a compor um grupo estatístico. Desse grupo estatístico, obtemos valores médios das grandezas macroscópicas de interesse, juntamente com seu erro estatístico. Dessa forma, temos a solução de sistemas macroscópicos através de simulações das interações microscópicas. Por se tratar de um método estatístico, é necessário que as simulações (também chamadas de histórias) sejam realizadas por diversas vezes até que o erro estatístico (desvio padrão) atinja a precisão desejada. Ou seja, enquanto a solução de um problema por métodos numéricos convencionais consiste na solução de equações algébricas e diferenciais, um problema simulado pelo Método de Monte Carlo tem sua física embutida nas PDFs, que definem a probabilidade da ocorrência dos fenômenos microscópicos, dessa forma, a execução da simulação resume-se ao uso de técnicas de amostragem para se chegar à solução desejada. Ao se simular o transporte de radiação, a história de uma partícula é descrita como uma sucessão de passos aleatórios seguidos de interações. Dependendo da interação sofrida pela partícula sua direção de propagação pode mudar, parte de sua energia pode ser depositada no meio e partículas secundárias podem ser geradas. A história de uma partícula é simulada até o momento em que ela abandone a geometria do problema, atinja uma energia mínima (energia de corte) ou seja absorvida. A história simulada de uma partícula não corresponde necessariamente ao comportamento de nenhuma partícula real, no entanto, seu comportamento é regido por distribuições estatísticas baseadas em modelos físicos. Desse modo, para um número grande o bastante de histórias, o resultado macroscópico obtido por análise estatística converge
21 4.2 Método de Monte Carlo 15 para o resultado correspondente à solução determinística da equação de Boltzman (que descreve o transporte de radiação). Embora o Método de Monte Carlo seja de mais fácil implementação que a solução da equação de Boltzman para problemas de geometria complexa, devido a sua natureza aleatória, o resultado obtido sempre terá um erro estatístico associado. Para reduzir esse erro, é necessário realizar simulações com um maior número de histórias, e isso implica em um maior custo computacional Técnicas de amostragem A amostragem correta de uma PDF é fundamental para que a simulação realizada pelo Método de Monte Carlo seja fiel ao processo representado pela PDF. Para garantir o caráter randômico da amostragem, um código de Monte Carlo emprega o uso de números racionais aleatórios 2 ξ uniformemente distribuídos entre 0 e Método da Inversão ou Método Direto Seja f(x) a função que descreve a densidade de probabilidade de ocorrência de um determinado fenômeno em um intervalo [a : b], sua função de distribuição acumulativa (CDF) é definida como sendo: x F (x ) = f(x)dx (4.1) a Sendo f(x) normalizada no intervalo [a : b], temos que F (x) é uma função contínua, crescente, com F (a) = 0 e F (b) = 1 que possui função inversa. Dessa forma, se temos F (x) = ξ, também temos uma função inversa F 1 (ξ) = x, que define um modo direto de se relacionar as variáveis x e ξ. Para a variável ξ variando uniformemente no intervalo [0 : 1] a probabilidade de se obter um determinado valor de x é igual a f(x). Exemplos: 2 Na prática, usa-se números pseudoaleatórios. Esse tema será abordado na seção 4.2.2, página 19
22 4.2 Método de Monte Carlo 16 Distribuição uniforme Tomemos como exemplo inicial o caso de uma distribuição uniforme já normalizada no intervalo [a : b]: f(x) = 1 b a Obtemos a CDF através da integração de f(x) no intervalo [a : x ]: F (x ) = x Isolando x tempos a função inversa: a 1 b a dx = x a b a = ξ (4.2) x = a + ξ(b a) = F 1 (ξ) (4.3) Função exponencialmente decrescente Para o caso de uma função exponencialmente decrescente, temos a distribuição normalizada no intervalo [0 : [ e obtemos a função inversa pelo mesmo método: f(x) = Σ T e Σ T x (4.4) x F (x ) = Σ T e Σ T x dx = 1 e Σ T x = ξ (4.5) 0 x = ln (1 ξ) Σ T = F 1 (ξ) (4.6) Embora a ideia de se calcular a função inversa seja conceitualmente simples, PDFs mais complexas podem fazer com que seu cálculo seja de difícil execução Distribuições Discretas O método da inversão também pode ser aplicado para distribuições discretas. Suponha que a variável x pode assumir qualquer valor inteiro entre 1 e N, com as respectivas probabilidades p1, p2,..., pn normalizadas.
23 4.2 Método de Monte Carlo 17 Desse modo, sua função cumulativa é dada por: 0 x < 1 F (x) = x i=1 p x 1 <= x <= N (4.7) 1 x > N Novamente, para um número aleatório ξ, de forma que F (x i ) < ξ < F (x i+1 ), o valor sorteado é x i. Tomemos como exemplo o lançamento de um dado de 6 faces. Sendo x o valor da face, ele tem igual probabilidade de assumir os valores 1, 2, 3, 4, 5 e 6. Dessa forma, sua PDF é definida por f(x) = 1/6,para x = 1, 2,..., 6. A figura 4.5 representa f(x). Figura 4.5: PDF para o lançamento de um dado de 6 faces Sendo assim, sua função cumulativa é definida por F (x) = x i=1 1/6, representada pela figura 4.6. Para um número aleatório ξ, obtemos o valor referente de x. O método discreto é amplamente aplicado em leitura de tabelas, no entanto, por depender do resultado de somatórias, pode ser mais lento que o método direto.
24 4.2 Método de Monte Carlo 18 Figura 4.6: Função acumulativa para um lançamento de um dado de 6 faces Técnica da rejeição O método da rejeição é um dos métodos mais simples de se amostrar uma PDF. No entanto, embora seja simples, apresenta uma eficiência menor do que dos outros métodos. O método consiste de duas etapas distintas. Na primeira, calcula-se o valor de f(x) para um valor de x aleatório variando uniformemente no intervalo a ser amostrado da PDF. a). Logo, se f(x) é definida no intervalo [a : b], calculamos f(x) para x = a+ξ(b A segunda etapa consiste em verificar se esse valor sorteado de x será aceito ou rejeitado. Para isso, definimos uma função auxiliar constante do tipo g(x) = c, de modo que g(x) f(x), para qualquer valor de x no intervalo [a : b]. Sendo ξ 2 um outro número aleatório, aceitamos o valor amostrado para x se: ξ 2 f(x) g(x) (4.8) Caso esta condição não seja satisfeita, um novo valor de x deve ser amostrado e o processo se repete.
25 4.2 Método de Monte Carlo 19 A figura 4.7 ilustra a técnica da rejeição. Figura 4.7: Técnica da Rejeição Embora a função g(x) possa ser definida com qualquer valor constante, contanto que seja sempre maior do que f(x), quanto mais próximas as duas funções forem, menor será a taxa de rejeição da técnica. O valor ideal para g(x) é ser igual ao valor máximo assumido por f(x) no intervalo em que é definida. A eficiência E do método que emprega a técnica da rejeição é dada pela razão entre as áreas definidas sob as funções f(x) e g(x), tal que: E = b a f(x)dx (b a)g(x) = 1 (b a)g(x) (4.9) Gerador de números pseudoaleatórios Um bom gerador de números pseudoaleatórios (do inglês Pseudo Random Number Generator, PRNG) é fundamental para um código que realize simulações através do Método de Monte Carlo. Uma sequência de números pseudoaleatórios é composta por números que si-
26 4.2 Método de Monte Carlo 20 mulam a aleatoriedade, no entanto são gerados por um algoritmo e possuem uma sequência bem definida. Uma sequência de números pseudoaleatórios apresenta um bom resultado quando nela aplicado um teste de aleatoriedade, mas possui a vantagem de poder ser repetida quantas vezes for necessário. Já uma sequência de números realmente aleatórios, não pode ser repetida, pois não é gerada por um algoritmo. A possibilidade de se prever o comportamento de uma sequência de números pseudoaleatórios e repeti-la é especialmente importante em simulações que usam o método de Monte Carlo por uma questão de reprodutibilidade. Em uma simulação de Monte Carlo que use números aleatórios, cada vez que a mesma simulação for realizada, um resultado diferente será obtido devido ao fato dos números sorteados não serem os mesmos. Obviamente, para um número de histórias suficientemente grande o valor estatístico tenderá a converter para o mesmo resultado, no entanto, cada história individualmente será diferente. A reprodutibilidade das sequencias de números pseudoaleatórios é importante quando se deseja realizar uma análise mais detalhada de uma simulação. Pois é possível saber exatamente a sequência de eventos que compôs cada história, ou então, ao se alterar algum parâmetro, saber o efeito que tal alteração trouxe na simulação, uma vez que a sequência de números deixa de ser uma variável do problema a ser analisada. Um dos métodos de se gerar números pseudoaleatórios é através de Geradores Lineares Congruenciais (GLC). O método consiste em um algoritmo que recebe uma semente (número inteiro positivo) como dado de entrada e retorna um número real entre 0 e 1 e uma segunda semente. Assim, o sistema pode ser retroalimentado e gerar uma sequência de números entre 0 e 1 que serão os números pseudoaleatórios. A estrutura de um GLC é representada pelo método recursivo, S n = (a.s n 1 + c)modm, (4.10) onde a, c e m são parâmetros fixos, S é a semente e a função (x)modm representa o resto da divisão de x por m. Depois de obtida a semente, o número pseudoaleatório é definido por, R = S/m (4.11)
27 4.2 Método de Monte Carlo 21 Vejamos um exemplo para o caso em que a = 47, c = 1, S0 = 1, m = 100 na tabela 4.1: Tabela 4.1: Exemplo de um gerador de números pseudoaleatórios de período curto S 0 = 1 S 1 = (47x01 + 1) MOD 100 = 48 MOD 100 = 48 R 1 = 48/100 = 0.48 S 2 = (47x48 + 1)MOD 100 = 2257 MOD 100 = 57 R 2 = 57/100 = 0.57 S 3 = (47x57 + 1)MOD 100 = 2680 MOD 100 = 80 R 3 = 80/100 = 0.80 S 4 = (47x80 + 1)MOD 100 = 3761 MOD 100 = 61 R 4 = 61/100 = 0.61 S 5 = (47x61 + 1)MOD 100 = 2868 MOD 100 = 68 R 5 = 68/100 = 0.68 S 6 = (47x68 + 1)MOD 100 = 3197 MOD 100 = 97 R 6 = 97/100 = 0.97 S 7 = (47x97 + 1)MOD 100 = 4560 MOD 100 = 60 R 7 = 60/100 = 0.60 S 8 = (47x60 + 1)MOD 100 = 2821 MOD 100 = 21 R 8 = 21/100 = 0.21 S 9 = (47x21 + 1)MOD 100 = 988 MOD 100 = 88 R 9 = 88/100 = 0.88 S 10 =(47x88 + 1)MOD 100 = 4137 MOD 100 = 37 R 10 = 37/100 = 0.37 S 11 =(47x37 + 1)MOD 100 = 1740 MOD 100 = 40 R 11 = 40/100 = 0.40 S 12 =(47x40 + 1)MOD 100 = 1881 MOD 100 = 81 R 12 = 81/100 = 0.81 S 13 =(47x81 + 1)MOD 100 = 3808 MOD 100 = 08 R 13 = 08/100 = 0.08 S 14 =(47x08 + 1)MOD 100 = 377 MOD 100 = 77 R 14 = 77/100 = 0.77 S 15 =(47x77 + 1)MOD 100 = 3620 MOD 100 = 20 R 15 = 20/100 = 0.20 S 16 =(47x20 + 1)MOD 100 = 941 MOD 100 = 41 R 16 = 41/100 = 0.41 S 17 =(47x41 + 1)MOD 100 = 1928 MOD 100 = 28 R 17 = 28/100 = 0.28 S 18 =(47x28 + 1)MOD 100 = 1317 MOD 100 = 17 R 18 = 17/100 = 0.17 S 19 =(47x17 + 1)MOD 100 = 800 MOD 100 = 0 R 19 = 00/100 = 0.00 S 20 =(47x00 + 1)MOD 100 = 1 MOD 100 = 1 R 20 = 01/100 = 0.01 S 21 =(47x01 + 1)MOD 100 = 48 MOD 100 = 48 R 21 = 48/100 = 0.48 Começando com a semente S 0 = 1, obtemos 48 como valor para a nova semente e 0.48 como primeiro número pseudoaleatório. Seguimos o processo até obtermos S 20 = S 0 = 1. À partir do aparecimento dessa semente pela segunda vez, a mesma sequência de números gerados passa a se repetir. Dizemos então que esse gerador possui um período de 20 números pseudoaleatórios. O período de um gerador depende, não só do algoritmo
28 4.2 Método de Monte Carlo 22 nele usado, mas também dos parâmetros que o definem. Se tivéssemos no caso anterior, S 0 = 4, a sequência seria S 0 = 4, S 1 = 89, S 2 = 84, S 3 = 49 e S 4 = 4 novamente, ou seja, um gerador com período de apenas 4 números. Para termos um gerador de longo período, seus parâmetros iniciais devem ser escolhidos de forma adequada. O primeiro parâmetro a se levar em consideração é o valor de m, módulo sob o qual as sementes são obtidas. Como as sementes são inteiras não negativas, o valor máximo de sementes possíveis e, consequentemente, números gerados é m 1. Sendo assim, geralmente usa-se o maior valor possível para m, que seria para sistemas computacionais de 32 bits e para sistemas de 64 bits RANECU Neste trabalho, o PRNG escolhido para ser usado nas simulações foi a adaptação para execução em GPU de um trabalho realizado por Badal em 2006 [4] que consistiu no desenvolvimento de uma versão paralela do RANECU com o intuito de permitir a execução em paralelo do código PENELOPE [4]. O gerador RANECU é usado em alguns softwares de Monte Carlo, dentre eles o PENELOPE, por possuir um período relativamente longo (2.3x10 18 números) e boas propriedades estatísticas. Tais características tornam o RANECU uma boa opção de gerador para códigos que não necessitam de sequências muito longas de números pseudoaleatórios. Desenvolvido por L Ecuyer em 1988 [19], combina as sequências produzidas por dois geradores congruenciais lineares multiplicativos (MLCGs) para produzir uma terceira sequência. Dadas as sequências R (1) i e R (2) i, a terceira sequência é definida por: R i = (R (1) i R (2) i )MOD(m (1) 1) (4.12) Cujo maior período possível é definido por (m (1) 1)(m (2) 1)/2. Para os parâmetros usados no RANECU (tabela 4.2), isso nos dá um período de números.
29 4.2 Método de Monte Carlo 23 Tabela 4.2: Parâmetros dos geradores MLCGs presentes no RANECU Módulo (m) Multiplicador (a) Gerador Gerador Existe também uma extensão do RANECU que trabalha com 3 geradores MLCGs para a qual o período é de aproximadamente 5x10 27 números. Nela temos m (3) = e a (3) = e a sequência final é definida por: R i = (R (1) i R (2) i + R (3) i )MOD(m (1) 1) (4.13) No entanto, por ser ligeiramente mais lenta, e as simulações de Monte Carlo tipicamente não precisarem de tantos números pseudoaleatórios, trabalharemos com a versão do código que possui apenas 2 geradores. código PENELOPE. A figura 4.8 mostra o código em FORTRAN do gerador RANECU presente no Paralelização do gerador de números pseudoaleatórios Os algoritmos para gerar números pseudoaleatórios baseados no método linear congruencial apresentam uma sequência de números com determinado período que passa a se repetir ciclicamente. Existem dois grandes grupos principais de geradores de números pseudoaleatórios em paralelo (PPRNG) que são os geradores que trabalham com diferentes sequências de números, independentes entre si, e aqueles que trabalham com diferentes trechos de uma única sequência numérica. Os geradores do primeiro tipo podem trabalhar com o mesmo algoritmo, mas com parâmetros iniciais distintos. Por sua vez, os geradores do segundo tipo podem ser divididos em dois subgrupos; os que tem um trecho da sequência destinada para cada gerador, e os que se utilizam da técnica conhecida como leapfrogging (saltos).
30 4.2 Método de Monte Carlo 24 Figura 4.8: Código em FORTRAN do gerador RANECU presente no código PENELOPE A técnica leapfrogging consiste em fazer com que cada gerador salte o número de geradores existentes entre um número e outro, ou seja, se possuo M geradores, cada gerador R i trabalhará com as sementes [S i, S i+m, S i+2m,...]. Desse modo, não existe problema de sobreposição entre os geradores. A outra técnica usada é a de segmentar a sequência original em diversas subsequências, deixando cada gerador responsável por um segmento. Se um PRNG possui uma sequência de N.M números pseudoaleatórios, ela pode ser separada em M geradores com N números cada. Assim, cada gerador R i trabalhará com as sementes [S in, S in+1, S in+2,..., S in+n 1 ]. Dessa forma, cada gerador poderá fazer uso de uma sequência única composta por N números. Naturalmente, caso esse valor seja excedido, o gerador invade a sequência do gerador seguinte. Deve-se, portanto, tomar cuidado ao dimensionar a paralelização de um PRNG para que a quantidade de números disponíveis
31 4.2 Método de Monte Carlo 25 por gerador seja o suficiente para realizar a simulação desejada. O PRNG presente neste trabalho é do tipo que trabalha com segmentos de uma sequência inicial. Para que seja possível seccionar a sequência de números gerada por um gerador em M subsequências, é necessário termos os valores das enésimas sementes [S 0, S n, S 2n, S 3n,..., S M.n ] pois elas servirão como parâmetro inicial para cada gerador, que trabalhará com uma sequência de n elementos. Da eq. 4.10, quando temos c=0, obtemos que o valor de S n é definido por: S n = (a n.s 0 )MODm, (4.14) Dessa forma, o valor de S n pode ser obtido sem a necessidade de se calcular o valor das sementes intermediárias, o que faz com que o processo de paralelização do gerador RANECU não seja difícil. Na seção 5.1 o método de paralelização do gerador RANECU é descrito em detalhes Transporte de fótons Na simulação do transporte de fótons, cada fóton é transportado individualmente. Sua história se inicia na fonte, ele realiza uma série de passos seguidos de interações, e tem sua história finalizada quando atinge uma energia mínima limite ou abandona a geometria do problema. Tipicamente, para energias iniciais da ordem de 1 MeV e meios com densidade próxima de 1g/cm 3, um fóton sofre até 10 interações antes de ter sua história finalizada. Desse modo, a simulação da história de um fóton é realizada passo a passo, ou seja, todas as interações e passos são simulados. Diferentemente do caso de elétrons, que sofrem centenas de interações, e portanto, normalmente cada interação simulada corresponde a uma média dessas interações. A simulação do transporte de um fóton é dividida nas seguintes etapas: Sorteio do passo 3 3 O passo de uma partícula é o caminho que ela percorre entre as interações
32 4.2 Método de Monte Carlo 26 Checagem das fronteiras Sorteio da interação Cálculo dos tallies 4 Essas etapas se repetem ciclicamente até que a história do fóton seja finalizada Sorteio do passo Ao se realizar o sorteio do passo de um fóton, deve-se sortear a direção de propagação, e o comprimento do passo. Usando coordenadas esféricas, isso significa sortear R, θ e ϕ em relação à posição atual da partícula, conforme ilustra a figura 4.9. Em seguida, esse valor é convertido para coordenadas cartesianas, seguindo a relação: S = ( x, y, z) = (R sin θ cos ϕ, R sin θ sin ϕ, R cos θ) (4.15) de modo que θ [0 : π] e ϕ [0 : 2π]. Para um sorteio isotrópico da direção de propagação das partículas, isto é, para que a projeção radial das partículas fique uniformemente distribuída sobre a superfície de uma esfera de raio R, a PDF utilizada para o sorteio da direção é dada por [25]: f(θ, ϕ)dθdϕ = 1 [ ] [ ] sin θ 1 sin θdθdϕ = 4π 2 dθ 2π dϕ (4.16) Isso nos mostra que θ e ϕ são variáveis independentes e podem ser amostradas pelo método da inversão. Desse modo, temos: θ = cos 1 (1 2ξ 1 ) cos (θ) = 1 2ξ 1 (4.17) ϕ = 2πξ 2 (4.18) 4 Os tallies são as informações macroscópicas calculadas a partir da simulação dos eventos microscópicos
33 4.2 Método de Monte Carlo 27 Figura 4.9: As coordenadas θ e ϕ descrevem a direção de propagação da partícula A figura 4.10 mostra a distribuição de pontos sob uma esfera de raio 1 seguindo a PDF anterior à esquerda e, à direita, a mesma distribuição caso o sorteio fosse feito segundo as funções θ = πξ 1 e ϕ = 2πξ 2. Podemos perceber que, embora seguir o sorteio de θ = πξ 1 pareça intuitivo, a distribuição de pontos decorrente deste sorteio apresenta um acúmulo nas regiões dos polos. Já para o sorteio de R, que descreve o comprimento do passo realizado pela partícula, a PDF que descreve o processo é f(x) = Σ T e Σ T x, onde Σ T é a seção de choque macroscópica do material no qual a partícula se propaga. Pelo método da inversão, obtemos: R = ln (1 ξ) Σ T = ln (ξ ) Σ T (4.19) A seção de choque macroscópica é uma grandeza que expressa o quanto um feixe de partículas é atenuado ao penetrar em um material. Para isso, ela relaciona as seções de choque microscópicas correspondentes a cada interação que a partícula pode
34 4.2 Método de Monte Carlo 28 Figura 4.10: Comparativo entre a distribuição dos pontos na superfície de uma esfera para o caso de uma distribuição isotrópica (à esquerda) e para um caso em que se amostra θ variando uniformemente entre θ e π (à direita) sofrer com o material. A figura 4.11 mostra as seções de choque de fótons para a água [6] Propagação da partícula A partícula pode ser propagada através de duas formas distintas. A primeira delas é o método convencional, onde a partícula é deslocada pelas diversas células da geometria, uma de cada vez. O segundo método é o delta tracking, mais conhecido como método de Woodcock, onde a partícula é propagada em um meio virtual homogêneo sem fronteiras. seguir. Ambos os métodos foram implementados nesse trabalho e serão discutidos à Método convencional (método com paradas) Pelo método convencional, após o passo da partícula ser sorteado, deve-se verificar se esse passo não implica em atravessar alguma das fronteiras que delimitam o meio no qual as partículas são propagadas. No caso da partícula cruzar alguma fronteira, ela é forçada a parar na interface e consultar as informações do novo material no qual ela está entrando.
35 4.2 Método de Monte Carlo 29 Figura 4.11: Seção de choque para fótons na água, considerando-se o efeito fotoelétrico, produção de pares, espalhamento coerente e incoerente Dessa forma, teremos três situações possíveis: Partícula não cruzou nenhuma fronteira Caso seja verificado que a partícula realizou o passo sem interceptar nenhuma fronteira, inicia-se a próxima etapa da simulação que é o sorteio da interação que ela irá sofrer. Partícula cruzou fronteira mas não mudou de material Caso seja verificado que, embora a partícula tenha interceptado alguma fronteira, o novo material continua sendo o mesmo (fronteira virtual), a partícula
36 4.2 Método de Monte Carlo 30 continua seu trajeto ignorando a existência dessa fronteira até que intercepte uma nova fronteira ou chegue à posição sorteada e, em seguida, sofra uma interação. Partícula cruzou a fronteira e mudou de material Nesse caso, a partícula terá um novo passo sorteado mas, dessa vez, usando as informações de seção de choque do novo material. No entanto, a direção de propagação não deve ser alterada, sorteia-se apenas o novo comprimento R do passo. Método de Woodcock ou delta tracking O Método de Woodcock [27] realiza a simulação do transporte da partícula ignorando as fronteiras existentes na geometria. O caminho livre médio (λ) de uma partícula depende do material no qual ela se propaga. Tal característica está relacionada à seção de choque macroscópica do material (Σ t ), de modo que, quanto maior o valor de Σ t menor o valor de λ. O método de Woodcock consiste em se considerar que a partícula é propagada em um meio homogêneo sem fronteiras composto pelo material do objeto simulador (meio) de maior Σ t, ou seja, aquele em que a partícula apresenta o menor caminho livre médio. Caso esse livre caminho médio seja maior do que a distância média entre as fronteiras, existe um ganho de desempenho. Para compensar a aproximação de meio homogêneo, é introduzida uma probabilidade de que a partícula sofra uma interação após realizar o passo (sem o método de Woodcock, considera-se 100% de chance de que haja a interação). A probabilidade de que ocorra uma interação é dada pela razão entre as seções de choque do meio (Σ Meio ) e a seção de choque máxima (Σ Máx ), ou seja; P rob = Σ Meio Σ Máx (4.20) Caso não haja interação, o passo é sorteado novamente até que ocorra. Método de Woodcock Vs método convencional Embora o método convencional obrigue a partícula a realizar uma parada toda vez que encontre uma fronteira, esse processo não é muito significativo (em termos de
37 4.2 Método de Monte Carlo 31 tempo computacional) para uma geometria composta por poucas fronteiras. Existe também a vantagem de sempre ocorrer uma interação após o passo ser concluído, e a existência de regiões com materiais de alta seção de choque não interfere na propagação da partícula em outras regiões da geometria. O método de Woodcock, por sua vez, mostra-se eficiente para geometrias compostas por muitas fronteiras, principalmente quando essas fronteiras não representam uma fronteira física real do meio. Para meios heterogêneos, o método de Woodcock ainda mostra-se eficiente quando as seções de choque entre os materiais não são muito diferentes. No caso de seções de choque muito diferentes, além das partículas serem propagadas com passos muito curtos, a probabilidade de interação também é baixa, fazendo que a eficiência do método caia Sorteio da interação Para o transporte de fótons, quatro interações principais são simuladas. O Efeito Fotoelétrico, Produção de Pares, Espalhamento Rayleigh (coerente) e Espalhamento Compton (incoerente), sendo que a probabilidade de ocorrência de cada uma dessas interações é proporcional à sua seção de choque. No decorrer das simulações de cada uma dessas interações, partículas secundárias como elétrons e pósitrons serão geradas. No entanto, o presente trabalho tem como proposta realizar apenas o transporte dos fótons. Desse modo, como aproximação inicial, sempre que uma partícula que não seja um fóton for gerada, sua energia será depositada localmente. Tal aproximação é justificada pelo fato do livre caminho médio do fóton ser ordens de grandeza maior do que o livre caminho médio do elétron. Dessa forma, mesmo que a simulação do transporte do elétron fosse realizada, haveria uma grande probabilidade de que o elétron acabasse depositando toda a sua energia antes de abandonar a célula (voxel) que foi gerado. Ainda assim, tal aproximação apresentará desvios de modelos físicos mais completos, por não considerar os fótons que serão gerados por esses elétrons. Fótons esses que provavelmente abandonariam o voxel e depositariam sua energia em outro local. Efeito Fotoelétrico
38 4.2 Método de Monte Carlo 32 Figura 4.12: Representação esquemática das 4 interações simuladas pelo código CUBMC No efeito fotoelétrico, o fóton interage com a eletrosfera do átomo sendo absorvido e transmitindo sua energia a um elétron que é ejetado. A conservação de momento se dá pelo recuo do átomo residual, e o elétron é ejetado com a energia do fóton (E) menos a energia de ligação que tinha com o átomo (U i ). Desse modo, para que seja possível a ejeção de elétrons mais fortemente ligados ao núcleo, é necessário que o fóton incidente possua maior energia. Quando um elétron mais fortemente ligado ao núcleo é ejetado, outro elétron mais energético transicionará entre as camadas eletrônicas para ocupar a vacância gerada. Tal transição gera fótons de fluorescência (ou raio-x característico) cuja energia é igual a diferença entre as energias dos níveis atômicos envolvidos nessa transição. Caso esse fóton seja capturado antes de abandonar o átomo, ele pode ejetar um segundo elétron chamado de elétron Auger. Nota-se então que a fluorescência e a emissão de elétron Auger são processos que competem entre si. A figura 4.13 contém uma representação esquemática do efeito fotoelétrico. A seção de choque fotoelétrica varia conforme a energia do fóton e o número atômico do átomo alvo. Ela é maior para fótons menos energéticos e átomos de maior número
39 4.2 Método de Monte Carlo 33 Figura 4.13: Representação do efeito fotoelétrico: O fóton incidente (1) excitando o átomo e causando a ejeção do elétron (3) gerando uma vacância (2). Em seguida ocorre o relaxamento atômico seguido da emissão de fluorescência (4) que pode, por sua vez, ejetar um segundo elétron (5), chamado de elétron Auger atômico. Conforme ilustra a figura Sendo um efeito predominante para fótons de baixa energia, a fluorescência decorrente do efeito fotoelétrico está, frequentemente, abaixo da energia de corte da simulação de Monte Carlo e não chega a ser simulada. Uma vez que o fóton de fluorescência é de baixa energia e o modelo físico simplificado usado nesse trabalho não engloba o transporte de elétrons, toda vez que a interação sorteada for o efeito fotoelétrico, toda a energia do fóton será depositada localmente e sua história será encerrada. Tal aproximação apresentará desvios significativos dependendo da energia do fóton. Os resultados obtidos serão apresentados e discutidos posteriormente. Produção de Pares A produção de pares é um fenômeno no qual o fóton é absorvido pelo átomo alvo dando origem a um par elétron-pósitron. Para que isso seja possível, o fóton incidente necessita ter no mínimo a energia de repouso do elétron e pósitron somadas,
40 4.2 Método de Monte Carlo 34 Figura 4.14: Seções de choque para o efeito fotoelétrico em função da energia para átomos de carbono, ferro e urânio 2m e c 2. pósitron. A energia restante é dividida entre as energias cinéticas do elétron e o E c = E γ 2m e c 2 = E c + + Ec (4.21) As equações que descrevem a seção de choque diferencial para a produção de pares são baseadas em modelos semi-empíricos. Neste trabalho é usada a tabela de seção de choque gerada pelo código PENELOPE que, para a produção de pares, é obtida por tabelas da XCOM [6]. Uma vez que as simulações deste trabalho não englobam transporte de elétrons ou pósitrons, toda vez que a interação selecionada for a produção de pares, toda a energia do fóton será depositada localmente e sua história será encerrada. As consequências de tal aproximação serão discutidas posteriormente. Espalhamento Rayleigh (coerente) O espalhamento coerente ocorre quando o fóton interage com a eletrosfera de um átomo mas não o excita ou ioniza. Ocorre a troca de momento entre o fóton e o
41 4.2 Método de Monte Carlo 35 átomo, mas ambos permanecem com a mesma energia. Tal espalhamento tem como análogo clássico uma colisão elástica entre bolas de bilhar. Para fótons incidentes com energia muito maior do que a energia de ionização da camada K do átomo, a seção de choque diferencial para o espalhamento Rayleigh por unidade de ângulo sólido pode ser aproximada por: dσ Ra dω = r2 e 2 (1 + cos2 θ)f 2 (q, Z) (4.22) onde θ é o ângulo de espalhamento, r e é o raio clássico do elétron e q é a magnitude da transferência de momento, definida por: q = p p = 2(E/c) sin(θ/2) (4.23) onde p e p são os momentos do fóton antes e após a interação e c é a velocidade da luz no vácuo. A função F (x, z) é o fator de forma atômico, que para átomos de densidade eletrônica esfericamente simétrica é dado por: F (q, Z) = 4π 0 ρ(r) sin ( ) qr r 2 dr (4.24) onde r é o raio, ρ(r) é a densidade em função do raio e é a Constante de Planck dividida por 2π. Descrevendo assim uma função monotonicamente decrescente em q, F (q, Z) varia de F (0, Z) = Z até F (, Z) = 0. Dessa forma, para fótons de baixa energia, a seção de choque para o espalhamento Rayleigh é dada por: qr E para fótons de alta energia: σ Ra 8 3 πr2 ez 2 (4.25) σ Ra E 2 (4.26) Neste trabalho, o algoritmo usado para a simulação do espalhamento coerente trabalha com as seções de choque obtidas pelo código PENELOPE.
42 4.2 Método de Monte Carlo 36 Espalhamento Compton (incoerente) O espalhamento Compton ocorre quando um fóton de energia E interage com um elétron que o absorve e emite um fóton secundário de energia E. Em seguida, o elétron é ejetado com energia cinética dada por E e = E E U i, onde U i é a energia de ligação do elétron. O análogo clássico do espalhamento Compton é o de uma colisão inelástica, onde o fóton transmite parte de sua energia para o elétron e é espalhado com um ângulo θ. No caso do espalhamento ocorrer com um elétron livre em repouso, para que haja conservação de energia e momento, é possível relacionar a energia final do fóton (E ) com o seu ângulo de espalhamento pela fórmula: onde k = E m ec. 2 E E 1 + k(1 cosθ) E c (4.27) A aproximação de Klein-Nishina para a seção de choque diferencial do espalhamento Compton para elétrons livres em repouso é dada por: dσ KN Co dω = r2 e 2 ( ) 2 ( Ec Ec E E + E ) sin 2 θ E c (4.28) As tabelas de seção de choque usadas nesse trabalho baseiam-se no modelo de Klein- Nishina, mas incluem fatores de correção para serem aplicáveis a elétrons ligados e em movimento. Assim como todas as tabelas usadas neste trabalho, essas também são as mesmas utilizadas pelo código PENELOPE [5]. O espalhamento Compton é o principal responsável por deposição de energia para fótons com o intervalo de energia de interesse desse trabalho (de 1 KeV à 9 MeV). Desse modo, o algoritmo que simula a interação Compton usado foi extraído do código PENELOPE sem simplificações no que diz respeito ao tratamento dos fótons. Os elétrons decorrentes do efeito Compton são gerados mas, em seguida, tem sua energia depositada localmente. Essa aproximação é necessária pois o código CUBMC não transporta elétrons. As consequências de tal aproximação serão discutidas posteriormente. Detalhes específicos da implementação do código serão apresentados na seção 5.
43 4.3 AMIGOBrachy Cálculo dos tallies Os tallies são os dados macroscópicos calculados através da simulação de eventos microscópicos da história de uma partícula, como fluxo de partículas, deposição de energia, densidade de partículas... Esse cálculo pode ser realizado para toda a geometria, ou apenas para algumas regiões de interesse. Nessa etapa de simulação da história do fóton, as informações pertinentes são transferidas para contadores e passam a compor o grupo estatístico que fará parte dos resultados da simulação. 4.3 AMIGOBrachy O software AMIGOBrachy foi desenvolvido para criar uma poderosa interface de fácil utilização que integra os planos de tratamento clínicos com simulações de Monte Carlo, fornecendo os principais recursos necessários para processar e editar imagens, importar e editar planos de tratamento e definir parâmetros para simulações em Monte Carlo analisando seus resultados. O modelo de implementação atual, usa o código MCNP6 para essas simulações, mas futuramente, também terá o CUBMC como uma opção de software simulador. O AMIGOBrachy foi desenvolvido usando a plataforma do MATLAB versão 8.0 (Mathworks Inc., Natick, MA), e tem como objetivo fornecer ferramentas eficazes para processar imagens médicas, identificar as posições de parada, detectar agulhas, criar objetos simuladores usando o formato DICOM, criar arquivos de entrada para realizar simulações de MC com o código MCNP6 e analisar os resultados através de uma interface de fácil utilização, que é composta por 8 módulos principais com ferramentas interativas específicas. O software é compatível com o sistema de planejamento de tratamento (TPS) Nucletron Oncentra T M (Elekta, Estocolmo, Suécia) e Varian BrachyVision T M (Varian Medical Systems, Palo Alto, CA), oferecendo a capacidade de importar planos de tratamento criados em ambos os TPS e realizar simulações utilizando o código MCNP6. As principais funcionalidades do software são: Importar ou criar uma sequência de imagens DICOM - atualmente imagens de tomo-
44 4.3 AMIGOBrachy 38 grafia computadorizada (TC) e ressonância magnética (RM) podem ser importadas, embora o usuário possa criar uma sequência de DICOM vazia e usar as ferramentas de desenho para criar objetos simuladores específicos, que podem ser exportados como imagens DICOM. Curvas automáticas de segmentação e de densidade de calibração não estão disponíveis para imagens de ressonância magnética, que devem ser segmentados usando contornos ou limites definidos pelo usuário; Elaborar e/ou importar estruturas - contornos criados com um TPS comercial podem ser importados e modificados, por exemplo, para excluir ou incluir uma região específica, e, em seguida, combinar com a imagem. Novos contornos também podem ser definidos; Importar ou criar um plano de tratamento - tempos de parada e posições podem ser editados ou criados pelo usuário ao visualizar uma distribuição de dose TG-43U1 antes de realizar uma simulação. Criar um arquivo de entrada para o MCNP6 definindo os parâmetros da simulação, como, a energia de corte, física detalhada ou simples, tallies e regiões de interesse. Além disso, pode-se escolher entre os seguintes opções de transporte: transportar os fótons na água e calcular a dose em água (D w,w ), transportar os fótons no meio e calcular a dose em água (D m,w ) ou transportar os fótons no meio e calcular a dose no meio (D m,m ). A energia média em cada voxel pode ser calculada. Esta parte do AMIGOBrachy pode ser facilmente adotada para preparar os arquivos de entrada para outros códigos de Monte Carlo; Visualizar e comparar os resultados, que são importados automaticamente no final da simulação e podem ser analisados utilizando o AMIGOBrachy ou exportados usando o formato DICOM que pode ser lido por TPS comerciais. Além das funcionalidades descritas, o AMIGOBrachy funciona como uma ferramenta auxiliar para criar arquivos de entrada à partir de geometrias do tipo mesh exportadas do Abaqus T M. O código CUBMC vem sendo desenvolvido visando uma integração futura com o AMIGOBrachy, como uma das opções de software para simulação do transporte de radiação.
45 5 Resultados 5.1 Gerador de Números Pseudoaleatórios (PRNG) O PRNG usado no código CUBMC é a versão adaptada para execução em paralelo do código RANECU. Tal procedimento foi explicado detalhadamente no artigo de Badal [4]. O algoritmo do gerador RANECU já foi apresentado na figura 4.8, pg.24, e a fórmula para se obter a n-ésima semente do gerador sem calcular as sementes intermediárias foi descrita pela eq. 4.14, pg.25. Sendo S i = (a i S 0 )MODm. Tal procedimento pode ser aplicado aos geradores MLCGs presentes no RA- NECU, pois ambos trabalham com o parâmetro c=0. O programa seedsmlcg foi desenvolvido para gerar essas sementes espaçadas usadas no código RANECU. Em seu artigo [4], Badal apresenta uma tabela (tabela 5.1) demonstrativa dos valores das sementes espaçadas, começando-se com S (1) 1 = S (2) 1 = 1 e promovendo um espaçamento de entre as sementes para diferentes geradores. O algoritmo implementado no código CUBMC foi baseado no código seedsmlcg e apresentou os mesmos resultados presentes na tabela 5.1 para as sementes S (1) i e S (2) i. As sementes S (3) i não são calculadas no código CUBMC. Após obter sementes concordantes com as da tabela 5.1, geradas pelo código seedsmlcg, o gerador de sementes presente no código CUBMC foi modificado para promover um espaçamento dinâmico entre as sementes, conforme a quantidade de geradores presentes em uma dada simulação. A sequência original de 2.3x10 18 números do gerador RANECU é dividida igualmente entre os geradores em execução no código. Desse modo, cada gerador trabalha com a maior sequência de números pseudoaleatórios possível, promovendo-se um aproveitamento máximo da sequência fornecida pelo gerador RANECU.
46 5.1 Gerador de Números Pseudoaleatórios (PRNG) 40 Tabela 5.1: Tabela demonstrativa do output do código seedsmlcg presente no artigo RANECU Seed 1 RANECU Seed 2 Extension Seed Comparativo de desempenho entre o RANECU e sua versão paralelizada Após a implementação da versão paralelizada do código RANECU, foi estabelecido um comparativo do tempo de execução do código em sua versão linear, com a sua versão para execução paralela. Para tal, a versão linear foi escrita em C, e é idêntica à versão presente no código PENELOPE mostrada na figura 4.8. A versão modificada para a GPU passa por algumas etapas extras. Inicialmente, as sementes para cada gerador são definidas por uma função executada fora da GPU, que é uma forma simplificada do programa seedsmlcg. Nessa função, foram definidas as sementes para 10 6 geradores afastadas entre si de Em seguida, essas sementes são transferidas para a GPU, onde os números aleatórios são gerados. Embora a execução em paralelo na GPU seja mais rápida, existe um tempo
47 5.1 Gerador de Números Pseudoaleatórios (PRNG) 41 computacional dispendido para preparar as sementes e transferir a memória que não existe na execução linear. Isso faz com que a execução em paralelo seja mais lenta para se gerar poucos números pseudoaleatórios. No entanto, a medida em que se aumenta a quantidade de números gerados, esse tempo inicial passa a se tornar desprezível frente ao tempo total, e o ganho de desempenho presente na execução em paralelo na GPU passa a ser evidente. Os números pseudoaleatórios foram gerados em um computador com Processador Intel Six-Core Core i7-980 de 3.33GHz e unidade de processamento gráfico NVIDIA GTX560Ti. A tabela 5.3 mostra o comparativo dos tempos entre as versões linear e em paralelo do gerador RANECU, onde o termo Ganho presente na tabela corresponde à razão entre esses tempos de execução. Tabela 5.2: Comparativo dos tempos entre as versões linear e em paralelo do gerador RANECU Números Gerados t Linear (s) t Paralelo (s) Ganho 1E x 1E x 1E x 1E x 1E x 1E x 1E x 1E x Podemos perceber que enquanto o método linear apresenta um tempo de execução diretamente proporcional à quantidade de números gerados, o método de execução paralela apresenta um tempo de latência de cerca de 0.15s. Tempo esse necessário para que as sementes sejam preparadas para os diversos geradores. Isso faz com que o método linear acabe sendo mais rápido para se gerar poucos números pseudoaleatórios. Após esse tempo de latência, a versão de execução paralela do RANECU gera cerca de 5.9x10 9 números por segundo, enquanto sua versão linear gera 5.8x10 7 números
48 5.1 Gerador de Números Pseudoaleatórios (PRNG) 42 por segundo, fazendo com que o ganho de desempenho ao se gerar uma maior quantidade de números pseudoaleatórios seja de cerca de 100x Função Rand Uma vez que o funcionamento adequado do gerador de sementes foi verificado, o gerador RANECU foi transferido para uma função do tipo device (que pode ser chamada durante a execução do kernel) chamada Rand. Essa função retorna um número pseudoaleatório entre 0 e 1 toda vez que for chamada. Os diferentes geradores são distribuídos entre os threads ativos na execução do código, de modo que cada thread gera uma sequência única de números. A figura 5.1 mostra o algoritmo correspondente à função Rand. Figura 5.1: Código em CUDA da função Rand Na figura, podemos observar que o código é muito semelhante ao código em
49 5.2 Cálculo de π 43 FORTRAN do gerador RANECU fig.4.8, pg.24. As principais diferenças são os termos idx, N, e o fato de que as sementes r.seed1[idx] e r.seed2[idx] são vetores na versão do código escrita em CUDA. O termo N designa o número de threads que serão executados no código. Esse número não significa necessariamente o número de threads ativos, mas sim o grau de paralelização da simulação. Para N thread s teremos N geradores diferentes. O termo idx representa o thread que está sendo executado no momento, variando de 0 a N. As sementes, por sua vez, são vetores pois cada thread acessará o valor de semente correspondente ao seu gerador. Mesmo que todos os geradores sejam executados simultaneamente, cada um deles acessará o valor de semente correspondente ao seu thread de índice idx. As variáveis I1, I2 e IZ são definidas no escopo de execução da função Rand, podendo ter seus valores apagados após a função ser executada. Isso permite que elas sejam alocadas na memória de registro dos threads, fazendo com que seu acesso seja mais rápido. No entanto, os valores das sementes devem ser preservados entre as chamadas da função Rand para que elas retornem o próximo número da sequência gerado pelo thread. Isso faz com que as sementes devam ser alocadas na memória global da GPU, de acesso mais lento. 5.2 Cálculo de π Tendo definido a função Rand, o passo seguinte é implementá-la em algum código simples para verificar se esta se comporta corretamente. Para isso, foi escolhido um algoritmo que calcula o valor de π pelo Método de Monte Carlo. Suponha uma circunferência de raio R = 1 inscrita em um quadrado de lado L = 2R = 2. A área da circunferência é dada por: e a área do quadrado é dada por: S o = πr 2 (5.1)
50 5.2 Cálculo de π 44 S q = (2r) 2 = 4r 2 (5.2) Desse modo, a razão entre a área do círculo e a área do quadrado é π/4. Nós não conseguimos determinar a área da circunferência de maneira exata sem saber o valor de π, mas se sortearmos um ponto aleatoriamente dentro do quadrado, sabemos que a probabilidade dele cair na área compreendida pela circunferência é proporcional à sua área. Dessa forma, se de um total de N pontos sorteados aleatoriamente, No caírem na região compreendida pela circunferência, temos que: N o N π/4 (5.3) Assim, o valor de π pode ser definido em função da razão entre os pontos que caem na área do círculo e o número total de pontos sorteados. Para que o algoritmo fique mais simples e, portanto, seja executado mais rapidamente, os pontos são sorteados em 1/4 da área total, apenas na região do primeiro quadrante da circunferência, conforme ilustra a figura 5.2. Figura 5.2: Representação gráfica de pontos distribuídos aleatoriamente para a determinação de π
51 5.2 Cálculo de π 45 Como R = 1, o sorteio de um ponto aleatório no quadrado resume-se a x = ξ 1, e y = ξ 2. Em seguida, temos dois contadores, um deles (N) que é acionado toda vez que um ponto é sorteado, e um segundo (N o ) que é acionado apenas se x 2 + y 2 < 1, ou seja, quando o ponto sorteado está na área do círculo. Finalmente, o valor de π é estimado por: π 4 N o N (5.4) Embora seja um problema bem simples, o cálculo de π pelo Método de Monte Carlo já incorpora os principais elementos presentes em simulações mais complexas. Além disso, como o valor de π é bem definido, é possível saber se a simulação está funcionando corretamente e avaliar a convergência de seu valor em função do número de pontos sorteados, mesmo após os números pseudoaleatórios começarem a se repetir. Também é possível fazer comparativos de análise de desempenho entre CPU e GPU. Em um primeiro teste, usamos 10 5 geradores em paralelo espaçados de sementes entre si, sorteando um total de 10 9 pontos. Dessa forma, cada gerador consumiu 2x10 4 números pseudoaleatórios (um para o sorteio de X e outro para Y ), estando longe da zona em que um gerador começa a invadir a sequência do outro. A tabela 5.3 apresenta a convergência para o valor de π em função da quantidade de pontos sorteados. Podemos notar que, embora apresente uma pequena oscilação, o valor da aproximação de π tende a convergir para o valor de π conforme o número de pontos sorteados aumenta. A figura 5.3 apresenta a convergência para o valor de π para os primeiros 1000 pontos sorteados em uma simulação contendo apenas um gerador pseudoaleatório. Em seguida, o mesmo teste de convergência foi feito para geradores com diferentes espaçamentos entre si. Foram usados 10 5 geradores para sortear 2x10 10 pontos, de modo que cada gerador ficou responsável por gerar 2x10 5 pontos, consumindo assim, 4x10 5 números pseudoaleatórios. O espaçamento entre os geradores foi dado por 10 K, onde k = 0, 2, 4, 6, 8, 10 e 12. Dessa forma, como cada gerador usará 4x10 5 números, para os valores de k = 0, 2 e 4, cada gerador começará a usar a sequência do gerador adjacente e é esperado que haja problemas de convergência.
52 5.2 Cálculo de π 46 Tabela 5.3: Convergência do valor de π Número de Pontos Aproximação de π Desvio % % % % % % % % % Figura 5.3: Convergência dos primeiros 1000 pontos para o valor de π em uma simulação com apenas um gerador pseudoaleatório conforme o esperado. A figura 5.4 mostra claros problemas de convergência para k = 0 e k = 2, A figura 5.5 apresenta uma vista aproximada das mesmas simulações, onde
53 5.2 Cálculo de π 47 Figura 5.4: Simulação de Monte Carlo para estimativa do valor de π usando diferentes configurações de geradores de números pseudoaleatórios. podemos observar que a simulação com k = 4 convergiu para um valor diferente de π. Enquanto as simulações com k > 5 (sem conflito entre os geradores) apresentam uma convergência adequada. Figura 5.5: Detalhe da simulação de Monte Carlo para estimativa do valor de Pi usando diferentes configurações de geradores de números pseudoaleatórios.
54 5.3 CUBMC 48 Conforme o esperado, as simulações com geradores de sementes iniciais pouco afastadas entre si apresentaram problemas de convergência ao se estimar o valor de π. Dessa forma, podemos concluir que a distribuição adequada das sementes entre os geradores apresenta uma grande influência no resultado das simulações. Por esse motivo, o código CUBMC adota um modelo dinâmico de espaçamento entre os geradores, a sequência original do RANECU é dividida igualmente entre o número de geradores presentes em uma dada simulação, garantindo assim, que cada gerador tenha a maior sequência de números possível antes de invadir a sequência do gerador adjacente. As diferenças nos tempos de execução das simulações em CPU e GPU seguiram o mesmo comportamento observado pelo gerador de números pseudoaleatórios. Para simulações com maior número de pontos sorteados (acima de 10 9 ) o ganho de desempenho da GPU foi próximo a 100x. 5.3 CUBMC Aspectos Gerais O código CUBMC (CUDA Based Monte Carlo) foi desenvolvido com o intuito de utilizar a capacidade de processamento paralelo das GPUs para atingir o maior ganho de desempenho possível em simulações de transporte de fótons através do Método de Monte Carlo. Para isso, o transporte é realizado inicialmente em geometrias de baixa complexidade, seguindo modelos físicos simplificados. Uma vez que o código apresente resultados satisfatórios dentro de suas limitações, seu grau de complexidade será aumentado a fim de obtermos um sistema mais compreensivo e versátil. Visando otimizações futuras, o código apresenta estrutura modular, permitindo que cada módulo possa ser substituído por um modelo mais complexo e abrangente em versões futuras.
55 5.3 CUBMC Geometria Desde seu primeiro modelo, o código CUBMC foi desenvolvido para simular o transporte de fótons em uma geometria descrita por voxels 5. Esse modo de descrever a geometria foi escolhido devido à sua simplicidade, mas que, no entanto, possibilita descrever facilmente geometrias complexas. A opção por voxels também facilitará a implementação de objetos simuladores provenientes de modelos anatômicos obtidos por imagens de tomografia computadorizada (CT), possibilitando assim, uma comunicação entre o código CUBMC e o sistema de planejamento AMIGOBrachy descrito na seção 4.3. No modelo atual do código CUBMC, os voxels são descritos por paralelepípedos que podem ter dimensões distintas em X,Y e Z, mas que, no entanto, devem ser todos do mesmo tamanho. O motivo dessa restrição está relacionado à rotina de checagem das fronteiras, e será abordado na seção No momento, existe essa restrição, mas será removida futuramente. Cada voxel da geometria é relacionado a um elemento inteiro de um vetor linear MAT que designa o material que o constitui. Dessa forma, é possível se determinar o material de cada voxel individualmente podendo, inclusive, cada voxel ser de um material diferente, dentro da limitação de memória da GPU para armazenar as tabelas de seção de choque desses materiais. O índice de cada voxel no vetor linear é dado por: ID = i + j DimX + k DimX DimY ; (5.5) onde i, j, k são, respectivamente os índices x, y e z do voxel, e DimX e DimY o número de voxels presentes em X e Y descrevendo a geometria. Como os índices dos voxels são iniciados em zero, em uma geometria contendo N voxels, estes serão identificados por ID variando de 0 a N-1, o que, tipicamente, implica em um vetor MAT contendo algo entre 10 6 e 10 9 elementos dependendo do problema. Apenas a memória global e a memória de textura são capazes de armazenar 5 Voxel é a extrapolação tridimensional do pixel. Do mesmo modo que um conjunto de quadrados (pixels) podem formar uma imagem bidimensional, um conjunto de paralelepípedos (voxels) definem uma geometria tridimensional
56 5.3 CUBMC 50 vetores com esse número de elementos 6. Atualmente, o vetor MAT encontra-se alocado na memória global, mas será movido para a memória de textura em futuras versões do código. Essa alteração deve resultar em uma leitura mais rápida do vetor, uma vez que a memória de textura possui cache de leitura dedicado. A versão inicial do vetor MAT contém apenas elementos com um valor inteiro que serve como índice para referenciar o material que constitui o voxel. Dessa forma, voxels contendo materiais de mesma composição atômica mas com diferentes densidades são referenciados como sendo materiais diferentes, e consultam tabelas de seção de choque distintas. A segunda versão do vetor MAT (que ainda será implementada) compreende materiais para os quais os valores das seções de choque macroscópicas variam linearmente com sua densidade, ou seja, para materiais de mesma composição atômica, o valor de suas seções de choque é diretamente proporcional à sua densidade. Assim, o vetor MAT será composto por elementos reais contendo, no mesmo número, informações sobre o índice do material e sua densidade. Dessa forma, materiais de mesma composição atômica consultarão a mesma tabela e terão, em seguida, os valores corrigidos segundo a sua densidade. Para um material de índice X e densidade ρ, o valor contido no vetor MAT será do tipo: MAT [ID] = 100X + ρ (5.6) pelas expressões: Desse modo, os valores do índice do material e sua densidade podem ser obtidos X = int(mat [ID]/100) (5.7) ρ = MAT [ID] 100X (5.8) onde int( ) é a função que retorna a parte inteira de um número racional. Naturalmente, é possível usar dois vetores distintos para armazenar as informações de material e densidade. No entanto, isso implicaria em um maior uso de 6 Os diferentes tipos de memória encontram-se descritos na seção , pg 10
57 5.3 CUBMC 51 memória da GPU, uma vez que seria necessário alocar um vetor de números inteiros (materiais) e outro de números reais (densidades), e também implicaria em duas leituras em vez de apenas uma. A GPU é significativamente mais lenta para realizar leituras do que operações, de modo que separar as informações de material e densidade de um único número é um processo mais rápido do que consultar duas informações alocadas em vetores distintos. Independentemente do formato do vetor MAT, ele não passa de uma forma de armazenar os índices referentes às tabelas de seção de choque a serem consultadas durante a execução do código Tabelas de Seção de choque As tabelas de seção de choque usadas no código CUBMC são geradas pelo código PENELOPE [5]. Entretanto, as tabelas geradas possuem mais dados do que são necessários para o código CUBMC, pois ele se restringe ao transporte de fótons, e apresenta diversas simplificações. Desse modo, da tabela originalmente gerada, o código carrega apenas o campo referente aos coeficientes de atenuação de massa para as interações simuladas (Espalhamento Coerente e Incoerente, Efeito Fotoelétrico e Produção de Pares) e alguns parâmetros a serem usados durante a rotina dessas interações 7. Os coeficientes de atenuação de massa (µ) tem unidade de cm 2 /g e representam a intensidade com que um fóton de uma determinada energia é absorvido por um material em decorrência de alguma interação. A tabela 5.4 apresenta o início da tabela com os coeficientes de atenuação de massa gerada para a água, pelo código PENELOPE e a figura 5.6 apresenta o gráfico da tabela completa. Na tabela 5.4, temos os coeficientes gerados para cada interação, e o coeficiente total que é a soma deles. Dessa forma, o coeficiente de atenuação de massa total representa a probabilidade de que o fóton sofra uma interação qualquer, e os coeficientes parciais representam a probabilidade de cada interação. O coeficiente de atenuação de massa total está relacionado ao livre caminho 7 Os parâmetros serão discutidos em detalhes quando o código que simula cada interação for abordado na seção 5.3.7
58 5.3 CUBMC 52 Tabela 5.4: Coeficientes de atenuação de massa gerada para a água, pelo código PENE- LOPE, para energias entre 50eV e 100eV. Onde 1mtu = 1g/cm 2 Energy Rayleigh Compton Photoelect. Pair Total (ev) (1/mtu) (1/mtu) (1/mtu) (1/mtu) (1/mtu) E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E E+04 médio da partícula e seu uso é descrito na seção Já o uso dos coeficientes parciais está relacionado ao sorteio de qual tipo de interação a partícula irá sofrer, e é descrito na seção Ao se executar o código CUBMC, os dados contidos nessas tabelas são transferidos para a memória constante da GPU, permitindo assim, um rápido acesso pelo thread. No entanto, a quantidade de memória constante presente na GPU é limitada a 64KB, de forma que, atualmente, o código CUBMC é capaz de simular uma geometria composta por até três materiais. Essa quantidade de materiais não é suficiente para simulações que tenham como objeto simulador um modelo anatômico proveniente de uma imagem de tomografia computadorizada. Dessa forma, em uma versão futura do código, será necessário armazenar as tabelas de materiais na memória de textura da GPU.
59 5.3 CUBMC 53 Figura 5.6: Coeficientes de atenuação de massa para a água gerados pelo código PENE- LOPE Modelos de fonte de radiação O módulo da fonte define onde os fótons são originados e qual a sua energia inicial, permite múltiplas fontes com formatos tridimensionais diferentes entre si. Também aceita que as energias sejam descritas por espectros, representando assim fontes reais. A descrição da fonte é realizada por um kernel à parte do kernel principal que realiza o transporte da radiação e as informações são transferidas por memória global. Dessa forma, antes que o kernel de transporte seja iniciado, a posição e energia inicial de cada partícula a ser transportada já foi definida pelo kernel da fonte. A fonte pode ser definida com qualquer forma desejada, pois não tem sua geometria vinculada à geometria do problema, podendo, inclusive, estar localizada pon-
60 5.3 CUBMC 54 tualmente na fronteira entre voxels. Distribuições de probabilidade espaciais, ou referentes ao espectro de energia também podem ser adotadas, no entanto, o modelo atual apresenta duas restrições. A primeira restrição é referente à localização espacial. O kernel de transporte finaliza a história de uma partícula quando ela se encontra fora da geometria descrita pelos voxels. Consequentemente, não é possível definir uma fonte externa à geometria do objeto simulador, pois a partícula já iniciaria sua história estando fora da geometria. Para se simular uma fonte externa, é necessário preencher o material entre a fonte e a região de interesse com vácuo. Tal processo é custoso, pois além de ocupar memória com informação de voxels que não são de interesse, as partículas serão transportadas nessa região. Esse processo é especialmente custoso devido ao fato de que todos os voxels da geometria são definidos por paralelepípedos de mesmo tamanho. O motivo dos voxels serem definidos assim será abordado na seção 5.3.6, mas o fato dos voxels serem todos do mesmo tamanho implica em que a partícula simulada leva quase o mesmo tempo computacional para ser transportada em um meio homogêneo de vácuo, quanto em uma geometria complexa. Essa forma de definir os voxels será alterada em versões futuras do código CUBMC mas, no momento, existe essa restrição. A segunda restrição referente à modelagem da fonte é que, por estar desassociada da geometria do objeto simulador, ela não inclui informações de material. Caso queira se definir algum material para a fonte, ele deve ser descrito juntamente aos outros materiais que compõem o objeto simulador e, portanto, ser definido por voxels. Embora a geometria da fonte possa ser descrita por qualquer função desejada, o material encontra-se vinculado ao modelo de voxels Sorteio do passo Após o kernel de transporte ser iniciado, a simulação da partícula passará ciclicamente por uma série de processos que descreverão sua história, até que ela abandone a geometria ou atinja a energia de corte. O sorteio do passo é o primeiro destes processos. O passo de uma partícula é o deslocamento que ela sofrerá entre duas interações
61 5.3 CUBMC 55 consecutivas e é realizado em coordenadas esféricas. Para isso, três valores devem ser definidos, R (Raio), θ (ângulo polar) e ϕ (ângulo azimutal). O sorteio de θ e ϕ é realizado segundo uma distribuição de probabilidade, que pode ser isotrópica, ou dependente da interação ocorrida, já o sorteio do Raio é realizado levando-se em consideração o livre caminho médio da partícula no meio e, para tal, depende do material no qual o fóton se propaga. A forma de se amostrar essas variáveis já foi demonstrada na seção pelas equações 4.17, 4.18 e Para que o Raio seja sorteado, é necessário consultar a tabela de materiais para se obter o valor da seção de choque macroscópica para fótons daquela energia, que é o coeficiente de atenuação de massa total multiplicado pela densidade do material. Embora a tabela gerada pelo PENELOPE forneça os coeficientes de atenuação de massa, ao ser transferida para a GPU, seu valor é multiplicado pela densidade do material para que o coeficiente de atenuação de massa total corresponda à seção de choque macroscópica. A figura 5.7 mostra o trecho do código que realiza a leitura da tabela e, em seguida, sorteia o valor do Raio. Figura 5.7: Trecho do código CUBMC que consulta a tabela de materiais e realiza o sorteio do comprimento do passo dado pelo fóton A primeira linha do código presente na figura 5.7 calcula o índice ID do voxel no qual a partícula se encontra. Em seguida, o loop for realiza a leitura da tabela do
62 5.3 CUBMC 56 material correspondente ao voxel, até obter qual é a linha (KE) da tabela, que é o índice do vetor no qual encotraremos as informações da tabela de seção de choque correspondentes à energia do fóton. Após o valor de KE ser determinado, os valores P 1, P 2, P 3 e P 4 correspondentes, respectivamente, às seções de choque parciais para o Espalhamento Rayleigh, Espalhamento Compton, Efeito Fotoelétrico e Produção de Pares são determinados, e a Seção de choque macroscópica (ST ) é dada como a soma desses termos 8. No caso da simulação ser realizada segundo o Método de Woodcock, é calculado o valor de ST para cada um dos materiais presentes na geometria e adotado o maior valor obtido. Após obter o valor de ST, a função RAND é executada para gerar o número pseudoaleatório r.rand[idx], e o valor do Raio é calculado conforme a fórmula Em seguida, é realizado o sorteio de cos θ e ϕ e a posição final (Xf, Yf, Zf) da partícula é determinada, conforme o código mostrado na figura 5.8. Figura 5.8: Trecho do código CUBMC que sorteia a posição (Xf, Yf, Zf) da partícula após o passo Da figura 5.8 podemos observar dois detalhes significantes deste trecho de código. O primeiro deles é que os ângulos θ e ϕ são definidos em relação a um eixo de coordenadas fixo. Como o sorteio da direção de propagação das partículas é feito na maioria dos casos de forma isotrópica, isso não é um problema. No entanto, após a partícula sofrer uma interação, é possível que exista uma nova distribuição de probabilidade para o 8 Embora a tabela gerada pelo código PENELOPE já forneça o valor da seção de choque total, é preferível economizar memória da GPU realizando a soma dos termos parciais
63 5.3 CUBMC 57 novo sentido de deslocamento, e essa distribuição trabalha com a referência no sentido de propagação que ela tinha antes de sofrer a interação. Nesses casos, é necessário realizar uma rotação no sistema de coordenadas. O segundo detalhe presente no trecho de código mostrado na figura 5.8 é o fato de existir um fator de escala para o Raio. Podemos observar que o valor do Raio é dividido por V x, V y e V z, que representam as dimensões em centímetros dos voxels em X, Y e Z, respectivamente. Esse fator de escala é necessário pois, durante a execução do kernel de transporte, os voxels são considerados sempre como sendo cubos de lado igual a 1cm. Isso faz com que a checagem de fronteiras seja mais simples, pois se o valor inteiro da posição da partícula mudou, significa que ela se encontra em outro voxel. No entanto, é necessário corrigir o valor do passo para que o deslocamento dos fótons no espaço dos voxels corresponda ao deslocamento no espaço real. Devido à essa aproximação usada para facilitar a checagem de fronteiras, todos os voxels devem ter as mesmas dimensões. Para que os voxels fossem de tamanhos diferentes, V x, V y e V z deveriam ser substituídos por vetores, correspondendo a diferentes fatores de escala para cada voxel. Essa alteração será implementada em futuras versões do código CUBMC Fronteiras Após realizar o passo, o segundo módulo do kernel de transporte verifica se a partícula não cruzou alguma fronteira. Caso a simulação esteja sendo realizada segundo o método de Woodcock, o módulo de checagem de fronteiras não é executado. Devido ao fator de escala usado que faz com que todos os voxels sejam tratados como cubos de lado 1cm, a condição que define se a partícula cruzou a fronteira de algum voxel é bastante simples. Basta verificar em X, Y e Z se a componente inteira da posição final e inicial é a mesma. Caso ela tenha mudado, significa que a partícula trocou de voxel. A figura 5.9 mostra o trecho de código que checa se houve colisão. A função floor(x) retorna apenas o valor inteiro de um número X. No entanto, não existe arredondamento, por exemplo, f loor( ) é 1, dessa forma, ao compararmos f loor(x) e floor(xf) não existe fatores de arredondamento por X ou Xf estarem próximos à borda
64 5.3 CUBMC 58 do voxel. Figura 5.9: Trecho do código CUBMC que checa se a partícula cruzou a fronteira do voxel ao realizar o passo Caso seja verificado que a partícula cruzou alguma fronteira ao realizar o passo, inicia-se uma rotina que calcula a posição X p, Y p e Z p, que correspondem às coordenadas em que a fronteira foi interceptada. A figura 5.10 representa o trecho do código que realiza a projeção de X p, Y p e Z p, para o caso da partícula estar situada em um voxel de índice X diferente. Figura 5.10: Trecho do código CUBMC que realiza a projeção de X p, Y p e Z p para uma partícula que se encontre em um voxel de índice X diferente do voxel de origem O processo de projeção realizado pelo código presente na figura 5.10 checa se a partícula avançou para um voxel de índice X maior ou menor, e corrige os valores para X p, Y p e Z p fazendo com que eles correspondam ao ponto de interceptação do voxel de índice X adjacente. Em seguida, o processo é realizado caso o voxel também possua índice
65 5.3 CUBMC 59 diferente em Y e em Z. Após o processo ser finalizado, os pontos X p, Y p e Z p correspondem ao ponto (na interface) em que a partícula cruza a fronteira do voxel inicial. A figura 5.11 representa esse processo de uma forma gráfica. Inicialmente, a partícula que estava na posição 1 tem seu passo sorteado para a posição 2. Após o primeiro ajuste ela tem sua posição projetada em 3, e em seguida, em 4, na fronteira do voxel. Figura 5.11: Representação gráfica bidimensional do processo de projeção para uma partícula que cruzou a fronteira entre voxels. A partícula que se encontrava na posição 1 teve seu passo sorteado para a posição 2. Em seguida, é projetada em 3 e em 4, sendo obrigada a parar na fronteira do voxel de origem. Feita a projeção na fronteira do voxel, o material entre os voxels adjacentes é consultado. Caso os voxels sejam do mesmo material, a partícula passa da posição 4 para a 3 e o material entre os voxels adjacentes é consultado novamente. Esse processo se repete
66 5.3 CUBMC 60 até que a partícula atinja a posição final, ou seja detectado que ela está entrando em um meio com outro material. Nesse caso, um novo valor de Raio deve ser sorteado usandose para o cálculo a seção de choque total do material no qual ela passa a se propagar. Apenas o valor do Raio é sorteado novamente, a partícula continua se propagando na mesma direção em que já vinha se propagando. Esse processo de checagem de fronteiras e novos sorteios de Raio é realizado até que a partícula atinja a posição sorteada. Nesse ponto deverá ocorrer uma interação Interações Após ser definido que a partícula sofrerá uma interação, é necessário sortear qual será a interação sofrida. Para isso, consultamos novamente os valores dos coeficientes presentes na tabela de seção de choque. Dessa vez, não será necessário ler as tabelas, pois o valor dos coeficiente já haviam sido lidos para se obter a seção de choque macroscópica, e ainda se encontram armazenados na memória de registro do thread. Conforme podemos observar no trecho de código presente na figura 5.12, ainda temos armazenados na memória o valor das variáveis ST, P 1, P 2, P 3 e P 4, que correspondem, respectivamente, à seção de choque macroscópica, e às seções de choque parciais para o Espalhamento Rayleigh, Espalhamento Compton, Efeito Fotoelétrico e Produção de Pares. Figura 5.12: Trecho do código CUBMC que realiza o sorteio entre as interações simuladas Para realizar a amostragem das interações, é utilizado o método discreto (descrito na seção ), onde o valor STS é uma variável pseudoaleatória uniformemente distribuída entre 0 e ST. Desse modo, a probabilidade de sorteio de cada uma das interações é proporcional à sua seção de choque. Dentre as 4 interações implementadas no código CUBMC, a única que não
67 5.3 CUBMC 61 apresenta simplificações (quando comparada ao código PENELOPE) é o Espalhamento Compton. A rotina que simula o Efeito Compton presente no código CUBMC é exatamente a mesma presente no código PENELOPE. A única alteração feita foi trocar os comandos goto presentes no PENELOPE por loops while. No entanto, a estrutura do algoritmo permanece a mesma. As outras interações foram inicialmente simplificadas, mas estão em processo de modificação para contemplar modelos físicos mais completos. No espalhamento Rayleigh, existe uma distribuição de probabilidade para a direção do fóton emergente. A aproximação atual adota uma distribuição isotrópica. No Efeito Fotoelétrico e na Produção de Pares, o modelo atual deposita toda a energia do fóton localmente e encerra sua história. Essa aproximação se justifica pelo livre caminho médio do fóton ser ordens de grandeza maior que o do elétron, fazendo com que a probabilidade do elétron abandonar o voxel seja baixa. No entanto, ignora a probabilidade de ocorrência de bremsstrahlung (efeito no qual um elétron é absorvido pelo material e ejeta um fóton), que é significativo em baixas energias. Essas aproximações serão corrigidas nas próximas versões do código CUBMC. No momento, o código tem seus testes realizados nos regimes em que o Efeito Compton é dominante Estimadores (Tallies) Os tallies correspondem às grandezas macroscópicas a serem obtidas como resultado da simulação. Basicamente, são contadores que armazenam as informações de interesse como fluxo ou dose nos voxels de interesse. Por se tratarem de dados obtidos estatisticamente, existe um erro estatístico intrínseco que também deve estar presente. Atualmente, existem dois tallies presentes no código CUBMC, o primeiro deles é o de cálculo de dose, que consiste em armazenar a informação a respeito da deposição de energia nos voxels, e o segundo tally é referente às interações ocorridas, quantas interações de cada tipo ocorreram e qual a energia por elas depositada. Este segundo tally é importante durante o desenvolvimento do código, pois permite estimar a participação de cada interação na deposição de energia total do sistema.
68 5.3 CUBMC Finalização da história da partícula A história da partícula é então descrita por uma sucessão de passos seguidos por interações até ser finalizada. Existem três condições que implicam na finalização da história da partícula. Ela abandonar a geometria do problema, atingir a energia de corte, ou atingir um valor limite de loops no contador. A figura 5.13 mostra o trecho de código no qual essas condições estão explícitas. Figura 5.13: Trecho do código CUBMC que define até quando a história de uma partícula é simulada As duas primeiras condições são estabelecidas devido ao seu significado físico. A história de uma partícula fora da geometria de interesse não é relevante para a simulação, e a energia de corte define o limiar inferior onde a energia depositada pela partícula é desprezível estatisticamente, ou os modelos físicos adotados não são mais válidos. O limite no contador de loops não possui significado físico, ele é apenas um processo preventivo para o caso de alguma partícula ficar em um loop infinito (muito mais demorado que os outros) e resultar em um atraso na simulação. O valor limite adotado é de 10 6, sendo que um fóton raramente atinge a marca de 100 loops antes de atingir sua energia de corte. O valor limite não interfere no resultado da simulação, apenas previne atrasos globais em decorrência de histórias pontuais de significância estatística desprezível Método de Woodcock O código CUBMC pode realizar a simulação da história da partícula usando o método convencional, que é o método de paradas, ou o Método de Woodcock, que foi descrito na seção A estrutura geral da simulação da história da partícula é praticamente a mesma, mas sem o módulo de colisão com as fronteiras. A partícula tem seu passo sorteado e se desloca ignorando qualquer fronteira do material.
69 5.3 CUBMC 63 Para o sorteio do passo, utiliza-se o valor da maior seção de choque macroscópica entre os materiais, que pode ser obtida de duas formas. A primeira delas consiste em ler a tabela de todos os materiais durante a execução do kernel, e escolher a de maior valor. A segunda forma seria realizar essa leitura fora do kernel, e gerar uma nova tabela contendo o valor de seção de choque a ser usada em cada energia. A vantagem do segundo método, é que, uma vez que a tabela seja gerada, a execução no kernel seria mais rápida, pois apenas uma tabela seria lida. No entanto, esse procedimento implicaria em se alocar mais uma tabela na memória constante da GPU, e isso não seria possível pois a memória constante tem apenas 64KB. Nas próximas versões do CUBMC, as tabelas de seção de choque serão alocadas na memória de textura e existirá uma tabela contendo apenas os valores de seção de choque máxima nas simulações usando o Método de Woodcock. No modelo atual, todas as tabelas são lidas e o maior valor é escolhido, conforme vemos no trecho de código presente na figura Com o valor de ST (seção de choque macroscópica) obtido, o passo é realizado ignorando a existência das fronteiras. No entanto, é introduzida uma probabilidade de que ocorra uma interação após o passo. Caso não ocorra interação, um novo passo é sorteado partindo da posição inicial. Esse processo é repetido até que ocorra uma interação. A probabilidade de ocorrência de uma interação é dada pela razão entre a seção de choque macroscópica do voxel no qual a partícula se encontra após o passo, e a seção de choque macroscópica máxima (usada para o sorteio do passo). Caso seja definido que ocorrerá uma interação, o processo se segue para a rotina das interações e tallies conforme o método de paradas Simulações As simulações de teste foram realizadas para quatro objetos simuladores distintos, compostos por água, osso, pulmão e osso+vácuo. Em cada um dos casos, a geometria do problema foi descrita por um objeto simulador cúbico de 64 x 64 x 64 cm 3, composto por 128 x 128 x 128 voxels, cada voxel com 0.5 x 0.5 x 0.5 cm 3. Foi simulada uma fonte de fóton pontual de 6 MeV posicionada em X = Y = Z = 32.1cm, contida no voxel de índice (64,64,64), posicionada ligeiramente deslocada do
70 5.3 CUBMC 64 Figura 5.14: Trecho do código CUBMC para se obter o maior valor de seção de choque macroscópica entre os diferentes materiais centro da geometria para não ser definida na fronteira entre voxels, pois tal posicionamento acarretaria em problemas para se realizar as mesmas simulações nos códigos de referência. As simulações foram realizadas através do Método de Paradas e de Woodcock e o valor da energia depositada nos voxels foi comparado com simulações realizadas com os códigos MCNP5 e PENELOPE. O tempo total de execução das simulações no código CUBMC foi comparado ao tempo de execução das mesmas simulações no código PENELOPE para se estimar o ganho de eficiência. Já para comparativos de valor de dose, foram usadas como referência as simulações realizadas com o código MCNP5. Essa escolha foi feita devido ao fato de termos acesso a um cluster para realizar as simulações em MCNP5. Dessa forma, foi possível simular um maior número de partículas de forma a minimizar o erro estatístico e obter
71 5.3 CUBMC 65 um melhor comparativo de dose. Embora não exista nenhuma fronteira física real em um material homogêneo, a geometria do problema é descrita por voxels, e essas fronteiras entre voxels obrigam as partículas a realizarem uma parada para consulta de material quando a simulação é realizada pelo método convencional. Isso gera um aumento significativo no tempo de execução, pois uma partícula pode cruzar centenas de fronteiras antes de ter sua história encerrada. As simulações realizadas usando o código PENELOPE foram executadas em um processador Intel Core(TM) i7-3520m CPU 2.90GHz, enquanto que as simulações com o MCNP5 9 foram executadas em um cluster SGI Altix XE 340 (composto por 96 processadores Intel Xeon six-core GHz), e as simulações do código CUBMC foram executadas em uma GPU GTX560Ti, da Nvidia. Para comparação de dose, foi escolhida uma coluna de voxels passando acima da fonte contida no plano médio entre a fonte e a fronteira da geometria. A coluna foi definida pelos voxels de índice V x = [0 : 127], V y = 64 e V z = 96, conforme ilustra a figura 5.15 Para um comparativo entre os tempos de execução, foram simuladas 2 x 10 8 partículas nos códigos PENELOPE e CUBMC. A tabela 5.5 apresenta os tempos de execução envolvidos nas simulações e o ganho de desempenho obtido pelas simulações em CUBMC frente as simulações realizadas no código PENELOPE. Tabela 5.5: Comparativo entre os tempos de execução para as simulações realizadas com o código PENELOPE e CUBMC PENELOPE CUBMC Woodcock Ganho Paradas Ganho ÁGUA 10051s 121s 83x 751s 13x PULMÃO 9839s 210s 47x 684s 14x OSSO 36030s 478s 75x 2349s 15x Podemos observar um significativo ganho de desempenho. No entanto, é necessário verificar se o valor de dose calculado está correto. Este comparativo não pode 9 Com exceção do caso da água, que foi simulado em um processador Intel Xeon E5645
72 5.3 CUBMC 66 Figura 5.15: Modelo fora de escala da posição da coluna de voxels (em azul) passando acima do voxel da fonte (em vermelho) ser feito de forma adequada para uma simulação contendo apenas 2 x 10 8 partículas, pois a flutuação estatística ainda é muito alta. Para um maior número de partículas, as simulações no PENELOPE levariam muito tempo (semanas). Desse modo, o código MCNP5 foi usado para comparativo de dose. A figura 5.16 mostra um comparativo de dose para as simulações em objeto simulador de água realizadas pelo CUBMC pelo método de paradas e de Woodcock, com a simulação realizada pelo MCNP5. Nas simulações realizadas através do código CUBMC, foram transportadas partículas, enquanto no MCNP5 foram simuladas 4 x Dessa forma, foi possível obter incertezas estatísticas inferiores a 1.5% nas simulações de ambos os códigos. A simulação em MCNP5 levou 18h31min para ser realizada em um processador Intel Xeon E5645 de 24 núcleos. Já as simulações com o CUBMC levaram 1h15min com o método de Woodcock e 7h38min com o método de paradas.
DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE ESTRUTURAS DE VOXELS USANDO UNIDADES DE PROCESSAMENTO GRÁFICO.
INSTITUTO DE PESQUISAS ENERGÉTICAS E NUCLEARES AUTARQUIA ASSOCIADA À UNIVERSIDADE DE SÃO PAULO DESENVOLVIMENTO DE UM SOFTWARE DE MONTE CARLO PARA TRANSPORTE DE FÓTONS EM ESTRUTURAS DE VOXELS USANDO UNIDADES