PARALELIZAÇÃO DO CÓDIGO DE MECÂNICA DOS FLUIDOS BRU3D
Leonardo Ostan Bitencourt IC
[email protected]
João Luiz Filgueiras de Azevedo PQ
[email protected]
Resumo
Este artigo tem por finalidade paralelizar o código BrU3D. Este código é capaz de simular
escoamentos de mecânica dos fluidos em torno de configurações aeroespaciais complexas de interesse do
Instituto de Aeronáutica e Espaço (IAE), como o Veículo Lançador de Satélites brasileiro (VLS). Quanto
maior a complexidade das geometrias simuladas, mais recursos computacionais como memória e,
principalmente, tempo de processamento são necessários. A paralelização do código BrU3D possibilitará
suprir estas duas limitações e representa um avanço no desenvolvimento das ferramentas de CFD do
Instituto de Aeronáutica e Espaço. Para tal, foi necessário um estudo minucioso de todo o processamento do
código BrU3D original para identificar todos os processos paralelizáveis. Como resultado, a solução do
código BrU3D paralelizado é idêntica ao do código serial, o que é mostrado nos resultados.
Abstract
The present paper has the objective of describing the strategy adopted for parallelization of the BrU3D
CFD code. This code is able to simulate aerodynamic flows over complex aerospace configurations of
interest to Instituto de Aeronáutica e Espaço (IAE), with particular emphasis on the configuration of the first
Brazilian Satellite Launcher (VLS). Complex configurations usually require additional computational
resources, namely memory and CPU time. The parallelization of the BrU3D code has the objective of
lessening such limitations by allowing the use of distributed memory PC clusters for the simulations of
interest. In this context, the work here reported represents an important advancement in the development of
CF tools at IAE. A careful analysis of the flow of information in the original BrU3D code was performed,
aiming at identifying the parallelizable processes. As a result, the solutions of the parallelized BrU3D code
are identical to those obtained with the serial code, as the simulations here reported indicate.
1.
INTRODUÇÃO
O campo da mecânica computacional dos fluidos (CFD)[1] já causou um grande impacto na ciência e
engenharia de dinâmica dos fluidos, com aplicações amplas que vão desde o projeto externo completo de um
avião até o projeto de motores mais eficientes. Contudo problemas de mecânica dos fluidos apresentam
dificuldade em sua resolução, mesmo em casos simples. A quantidade de problemas de engenharia que não
apresenta solução analítica supera em muito os que as têm. Por isso, antigamente, engenharia se fazia
principalmente por exaustivos ensaios em laboratório. Após o surgimento do computador, programas
começaram a auxiliar o engenheiro em todo o seu projeto, desde a concepção estrutural em ferramentas CAD
até a fase de testes.
O trabalho desenvolvido se insere no desenvolvimento de ferramentas numéricas de mecânica dos
fluidos computacional que simulem escoamentos em configurações aeroespaciais de interesse para o Instituto
de Aeronáutica e Espaço (IAE/CTA), sendo desenvolvido pelo grupo de CFD deste instituto. A meta do
grupo de CFD do IAE é o desenvolvimento de ferramentas computacionais utilizadas no estudo das
propriedades aerodinâmicas e aerotermodinâmicas de problemas de interesse aeroespacial, como veículos de
sondagem e do primeiro veículo lançador de satélites brasileiro (VLS).
O código BrU3D[2-3] resolve numericamente as equações de Navier-Stokes para o regime estacionário
através do método de volumes finitos em malhas não-estruturadas. Ele apresenta precisão de 2ª. ordem no
tempo e espaço, pois utiliza um método de Runge-Kutta de 5 estágios adaptado para avançar no tempo. A
implementação foi feita na linguagem de programação Fortran 77 e apresenta quase 22500 linhas de código.
O objetivo deste trabalho foi implementar uma solução eficiente e barata que aumentasse a capacidade
de processamento do laboratório computacional do IAE, expandindo a memória disponível e reduzindo o
tempo necessário para a resolução dos problemas de CFD. A solução encontrada foi utilizar
microcomputadores já existentes ligados em rede e distribuir o processamento que antes só se dava em uma
máquina. Para distribuir o processamento entre os computadores [4], ou nós, utilizou-se da biblioteca MPI[5-7]
(Message Passing Interface), além da biblioteca METIS para particionar de maneira muito eficiente grafos
não-estruturados.
2. ASPECTOS TEÓRICOS E NUMÉRICOS
A Mecânica Computacional dos Fluidos (CFD) pretende desenvolver a capacidade de simular o
comportamento de um fluido tridimensional, compressível, viscoso e turbulento. Para isso, é necessário
modelar matematicamente e numericamente o fluido através das equações de Navier-Stokes, que serão
descritas abaixo.
2.1 FORMULAÇÃO TEÓRICA DE MECÂNICA DOS FLUIDOS
Estas equações representam a formulação mais geral para um escoamento de um fluido contínuo. As
equações de Navier-Stokes, para um gás perfeito, sem a geração de calor e com forças de campo
negligenciáveis podem ser escritas conforme apresentado na Eq. 1-3.
dfsdfsdfdsfdsfsdfsdfsdfsdfsdfsdfsdfsdfsadfasdfdsfsdfd(1)
dfsdfsdfdsfdsfsdfsdfsdfsdfsdfsdfsdfsdfsadfasdfdsfsdfd
sdfdsdsdfsdfsdfsdfasdfsdfsdfsadfsadfsadfsadfsdafasdf(2)
vsdfsdfsdfsdfdsfsdfsdfsfsadfsads fsadfsdfsdfsdfsdfsadf(3)
Nas equações acima, ρ, p e H são, respectivamente, a densidade, a pressão e a entalpia total ou de estagnação;
é a velocidade do fluido, é o tensor viscoso, é o vetor do fluxo de tranferência de calor e t é o tempo.
Este sistema de equações não está fechado, ainda são necessárias mais equações para que o número de
variáveis seja igual ao de incógnitas. No contexto do código BrU3D, o fluido pode ser aproximado a um gás
perfeito. Então, a equação para um gás perfeito pode ser adicionada ao sistema,
(4)
p = ρ R T = (γ – 1) ρ ei,
onde a enegia interna específica, é dada por
ei = Cv T .
(5)
Aqui, Cv é o calor específico à volume constante, T é a temperatura do fluido e γ = Cp/Cv é a razão dos
calores específicos. O vetor de fluxo de calor é calculado de acordo com a lei de Fourier para a condução de
calor como
,
(6)
onde k é o coeficiente de condutividade térmica do fluido. As componentes do tensor viscoso são dadas por
ddddddddddddddddddddddddddddddddddd(7)
onde µ é o coeficiente de viscosidade do fluido e λ é o segundo coeficiente de viscosidade. Para as aplicações
de interesse, a hipótese de Stokes pode ser assumida, que fornece
λ = -2/3 µ .
(8)
A entalpia total ou de estagnação de um fluido é definida como
2
,
(9)
e a derivada substancial é dada por
.
(10)
2.2 COMPUTAÇÃO PARALELA
O processamento paralelo[4] permite a solução de problemas que exigem muito processamento e
memória utilizando computadores de modestos recursos computacionais, pois decompõe processamento do
problema entre vários computadores. Assim, contribui para o desenvolvimento da ciência, seja por reduzir
limitações de memória, seja por aumentar a velocidade de processamento, ou na maioria dos casos ambas.
No presente trabalho, a arquitetura distribuída foi adota, pois era necessário expandir a capacidade de
processamento a um custo baixo. As arquiteturas de memória compartilhada apresentam a grande vantagem
de serem mais rápidas, pois não há necessidade da troca de dados e o acesso ao bus é centenas de vezes mais
veloz. Contudo, o custo de cada supercomputador é altíssimo, sendo da ordem de alguns milhões de dólares.
A decomposição do problema, independente da forma como é feita, deve priorizar um balanço
adequado de carga e a minimização de comunicação entre processadores no caso de memória distribuída. O
balanço de carga adequado consiste em dividir igualmente o trabalho entre os processadores disponíveis. Já a
minimização da comunicação entre os nós, como o nome já indica, visa diminuir ao máximo a comunicação
ou troca de mensagens, pois representam tempo a menos de processamento que seria gasto na solução do
problema.
2.3 PARTICIONAMENTO DE MALHAS COMPUTACIONAIS
A malha computacional representa um espaço discretizado descrito por células tridimensionais
formadas por pontos em coordenadas cartesianas. A paralelização do código BrU3D utiliza o conceito de
decomposição de domínios, ou seja, cada nó do cluster resolve as equações de Navier-Stokes para somente
algumas partições de células. Por isso, é necessário previamente dividir/particionar o domínio/malha
computacional para que todos os processadores atuem em seus respectivos domínios. Para isso, utilizou-se
uma biblioteca robusta de particionamento de malhas e grafos chamada METIS. O resultado do
particionamento de uma malha não-estruturada pode ser visualizado na Fig. 1.
Figura 1 : Visualização de uma malha computacional particionada pelo METIS.
3.
SOLUÇÃO IMPLEMENTADA NO CÓDIGO BRU3D
Após estudar detalhadamente o funcionamento do código BrU3D, e pesquisar na literatura trabalhos
semelhantes de paralelização, chegou-se à conclusão que a melhor estratégia de paralelização seria efetuar a
decomposição de domínios ou decomposição da malha computacional. Adotando esta estratégia, alterações
mínimas no código BrU3D original foram necessárias, pois cada partição ou domínio decomposto é tratado
como uma malha computacional. Nas interfaces entre domínios é criada uma nova condição de contorno, que
é atualizada com as propriedades das partições vizinhas. Logo, o projeto possui duas etapas:
3
•
•
Decomposição da malha computacional;
Gerenciamento da comunicação entre os processadores.
3.1.
DECOMPOSIÇÃO DA MALHA COMPUTACIONAL
A decomposição da malha, ou domínio computacional, consiste em gerar N subdomínios que, juntos,
reconstituem o domínio original. Uma decomposição bem feita, considerando-se o mesmo tempo de
processamento por volume, deve gerar subdomínios com aproximadamente o mesmo número de volumes. O
resultado é um array de tamanho igual ao número de vértices do grafo que carrega como informação a relação
nó-partição. Com esta informação, dois arquivos são gerados. Um contendo a volume-partição e outro
contendo as faces de interface entre cada partição.
3.2.
COMUNICAÇÃO ENTRE OS PROCESSADORES
A condição de contorno de interface representa uma extensão do subdomínio. Esta condição de
contorno é implementada através de volumes ghost atualizados, a cada iteração, com as propriedades do
volume pertencente ao subdomínio vizinho e que compartilha a face onde foi aplicada a condição de contorno
de interface. A Fig. 2 ilustra o funcionamento da condição de contorno.
Figura 2: O volume ghost p na partição 1 possui as propriedades do volume p na partição 2, já o
volume ghost k na partição 2 possui as propriedades do volume k da partição 1.
4. VISÃO GERAL DO FUNCIONAMENTO DO CÓDIGO BRU3DPAR
Cada processador ou nó do cluster possui uma cópia do programa BrU3DPar, assim como uma cópia
dos arquivos de malha.. Após o comando “mpirun”, o programa BrU3DPar é acionado simultaneamente em
todos os computadores/nós do cluster. A inicialização do código Bru3dPar diferencia-se muito pouco da
inicialização do código BrU3D original. Os arquivos de malha anteriores são lidos normalmente, contudo,
cada processador somente armazena os volumes e condições de contorno do seu domínio. Isso é possível
através da leitura do arquivos que contém as relações volume-partição e as condições de contorno de
interface. A cada iteração do código, o resíduo máximo e as condições de contorno são autalizadas, inclusive,
as novas condições de contorno de interface.
O pseudo-algoritmo abaixo mostra o funcionamento do BrU3DPar.
Início
Inicialização do MPI ;
Carrega malha computacioanl;
Atualiza Condições de Contorno de Interface; // INTERFACE
Enquanto (A solução não Converge OU número de iterações não excedeu o limite)
Início
Avança no tempo; // Aplicação do Runge-Kutta de 5 passos
Atualiza condições de contorno;
Se (Necessário SALVAR a solução atual) Então
Início
Gera o arquvio TECPLOT
Gera o arquivo com a solução atual
Fim
Fim
4
5.
RESULTADOS
O resultado apresentado neste relatório foi simulado num cluster formado por dois nós ligados a uma
rede 100Mbits. A solução numérica convergida apresenta máximo resíduo da densidade no escoamento da
ordem de 10-7 . Compara-se o resíduo do código BrU3D original com o código BrU3DPar, faz-se o mesmo
com a distribuição de Cp ao longo do hemisfério-cilindro. Também é exibida uma visão 2D da configuração
simulada. O resultado apresentado não tem por objetivo comparar desempenho, apenas mostrar a coerência
dos resultados numéricos.
A malha computacional mostrada na Fig. 3 é uma malha tridimensional sobre a configuração de um
hemisfério-cilindro. Ela é composta por 21916 tetraedros e é considerada uma malha modesta em relação às
malhas nas quais se pretende utilizar o código BrU3DPar. Especificamente, a Fig. 3 apresenta apenas a malha
superficial, composta evidentemente de triângulos.
Figura 3: Visão da malha superficial sobre o hemisfério-cilindro
Os parâmetros de simulação utilizados foram número de CFL=0.2, número de Mach do escoamento
não perturbado, M∞ = 2.0 e ângulo de ataque α = 0 graus.
O decaimento do resíduo máximo exibido na Fig. 4(a) demonstra que as duas soluções numéricas estão
convergindo, ou seja, a solução numérica do escoamento está a cada iteração se aproximando da solução
exata. Além disso, as duas curvas sobrepuseram-se exatamente mostrando que a paralelização não alterou em
nada o processo de cálculo do código original. O gráfico da Fig. 4(b) demonstra a coerência dos resultados.
Figura 4: (a) Gráfico comparativo do decaimento do resíduo máximo pelo número de iterações do
BrU3D Serial e do BrU3DPar Paralelo. (b) Gráfico comparativo de Cp ao longo da superfície do hemisfério.
5
Como é possível observar na Fig. 4(a), os resultados do BrU3D original se sobrepuseram aos valores
calculados pelo BrU3DPar, como não poderia ser diferente. Logo abaixo, na Fig. 5 é apresentada uma visão
2D do campo de velocidade e pressão em torno do hemisfério-cilindro, mostrando contornos de pressão e
vetores velocidade em um plano longitudinal do escoamento e sobre a superfície do corpo.
Figura 5: Campo de velocidades ao longo do hemisfério-cilindro além das curvas de pressão
6.
CONSIDERAÇÕES FINAIS
A paralelização do código BrU3D foi realizada com sucesso, os resultados numéricos foram
consistentes e o resíduo convergiu como esperado, mostrando a consistência da paralelização e das bibliotecas
utilizadas, como LAM-MPI e Metis. Contudo, devido à falta de tempo e máquinas suficientes para realizar
testes de desempenho, não foi possível medir o ganho de desempenho natural decorrente da paralelização. Na
próxima bolsa de iniciação científica, o trabalho de paralelização terá continuidade e avaliações de
desempenho serão feitas.
AGRADECIMENTOS
Os autores agradecem ao CNPq pelo suporte fornecido através do Processo No. 501200/2003-7 e do
Processo No. 502053/2004-6. O primeiro autor gostaria ainda de agradecer aos colegas Edson Basso e Enda
Bigarella pela ajuda fornecida, que se mostrou essencial para o bom desenvolvimento do presente projeto.
REFERÊNCIAS BIBLIOGRÁFICAS
[1] Fletcher, C.A.J., Computational Techniques for Fluid Dynamics, Vol. I, Spring-Verlag, New York,
1988.
[2] Scalabrin, L.C., “Numerical Simulations over Three-Dimensional Flows over Aerospace
Configurations”, Master Thesis, Divisão de Engenharia Aeronáutica, Instituto Tecnológico de
Aeronáutica, São José dos Campos, SP, 2002.
[3] Bigarella, E.D.V, Basso, E., and Azevedo, J.L.F., “Multigrid Adaptative-Mesh Turbulent Simulations
of Launch Vehicle Flows”, AIAA Paper No.2003-4076, Proceedings of the 21st AIAA Applied
Aerodynamics Conference, Orlando, FL, June 2003.
[4] http://www.cs.uncc.edu/~abw/parallel/par_prog/ - “Parallel Programming: Techniques and
Applications” – UNC Charlotte.
[5] http://www.ncsa.uiuc.edu/Divisions/eot/Training/ - Tutorial on-line de MPI do National Center for
Supercomputing Applications (NCSA).
[6] http://www.lam-mpi.org – Site oficial da biblioteca LAM-MPI..
[7] http://www-unix.mcs.anl.gov/mpi/learning.html - Tutorial on-line de MPI.
6