Simulação eficiente de fluidos no espaço paramétrico de malhas estruturadas tridimensionais

1. Tese de Doutorado Simulação eficiente de fluidos no espaço paramétrico de malhas estruturadas tridimensionais Aluno: Vitor Barata R. B. Barroso (vbarroso@inf.puc-rio.br) Orientador: WaldemarCeles(celes@inf.puc-rio.br)

2. Motivação • Exemplos de fluidos • Correntes de rios • Ondas do mar • Fumaça de uma chaminé ou cigarro • Vapor quente de um bule ou motor • Estudos importantes envolvendo fluidos • Vento passando pelas asas e turbinas de um avião • Formação e evolução de tempestades e furacões • Fluxo de sangue por veias e artérias • Transporte de água e óleo em dutos

3. Introdução • Simulação de fluidos • Métodos Lagrangianos: representam o fluido com partículas (ex.: SPH) • Métodos Eulerianos: subdividem o espaço em voxels • Descrição euleriana de um fluido em uma grade regular • Campos de valores escalares e vetoriais amostrados no centro de cada célula (i,j,k) da grade • Posição: x = [x y z]T • Velocidade: u = [u v w]T • Pressão: q • Propriedades intrínsecas constantes • Massa específica:  • Viscosidade cinemática:  • Temperatura: T • Objetivo: integrar os campos de propriedades para descrever sua variação no tempo x(i,j,k) y x

4. Fronteiras Curvas • Grade regular • Marcação de células exteriores • Tratamento diferenciado de células de fronteira • Grade em multirresolução • Captura melhor a forma • Não evita artefatos associados a degraus • Malha simplicial • Requer informação topológica • Maior custo computacional • Grade curvilínea paramétrica estruturada • Simples e rápido • Não requer topologia • Evita artefatos

5. Proposta • Idéia básica: explorar o conceito de malhas curvilíneas de CFD • Mapear a grade física curvilínea do espaço do mundo para uma grade computacional uniforme no espaço paramétrico • Realizar a integração eficientemente no espaço paramétrico • Vantagens: • Evitar a necessidade de armazenar dados para toda a caixa envolvente do domínio • Considerar paredes e obstáculos curvos de maneira natural, sem a necessidade de se refinar a grade ou de se lidar com casos especiais próximo às fronteiras

6. Proposta • Foco: • Gerar boas aproximações para aplicações interativas e de visualização científica • Estratégia: • Começar com o algoritmo Stable Fluids clássico [Stam99] • Realizar adaptações para levar em conta métricas e Jacobianos, além do esquema de amostragem deslocada • Incorporar gradativamente técnicas avançadas mais recentes conforme o necessário, sempre verificando as adaptações necessárias para a sua implementação no espaço paramétrico

7. Transformação de Coordenadas • Regra da Cadeia para gradientes de campos escalares • Notação Matricial • Jacobiano

8. Transformação de Coordenadas • Regra da Cadeia para velocidades • Notação Matricial • Jacobiano

9. Transformação de Coordenadas Domínio físico(espaço do mundo) Domínio computacional(espaço paramétrico) • J e J-1 podem ser precomputados • armazenados por célula como 9 escalares (4 em 2D) • podem ser interpolados fora dos pontos amostrais • Precisaremos também de termos de segunda ordem: • sxx, syy, szz,txx, tyy, tzz,pxx, pyy, pzz

10. Modelo para Simulação • Equações de Navier-Stokes (fluidos incompressíveis) • Termos: • Advecção: as propriedades do fluido são carregadas por ele próprio • Pressão: o fluido tende a ocupar áreas de menor pressão • Difusão / viscosidade: atrito interno e resistência à deformação • Forças externas: acelerações adicionais

11. Algoritmo Stable Fluids original • Equações de Navier-Stokes (fluidos incompressíveis) • Stable Fluids [Stam99] • Aproximação de alto desempenho e incondicionalmente estável • Separação de primeira ordem em passos fracionados • Termos considerados sequencialmente de forma independente • Algoritmos eficientes e estáveis para cada etapa • Cálculo da velocidade em cada passo usa o resultado do anterior • Ordem dos operadores é importante • Advecção exige campo livre de divergência

12. Separação de Operadores • Estratégia do tipo “dividir-e-conquistar” • Dividir um problema complexo em um conjunto de subproblemas mais simples • Resolução de cada subproblema por algoritmos especializados • Soma das contribuições de cada solução intermediária • Tipos • Separação de coordenadas • Passos fracionados (separação por fenômenos físicos) • Erro • Pode ser estimado, existem métodos de diferentes ordens • Pode limitar a acurácia da integração mesmo com operadores exatos

13. Separação de Operadores • Problema de valor inicial: • Resolução direta: • Método de Lie-Trotter (primera ordem):

14. Discretização do domínio • Amostragem colocalizada[Stam99] • Pressões e velocidades nos centros das células • Equações desacopladas, componentes oscilatórias • Amostragem deslocada[Harlow65,Bridson08] • Pressões nos centros, velocidades nas faces das células • Componentes de velocidade separadas (interpolação)

15. Representação das velocidades Velocidades no mundo • Velocidades no espaço do mundo • Tratamento matemático mais simples • Malha deslocada: • Desalinhamento entre velocidades e normais das faces não é prejudicial ao resultado[Azevedo12] • Maior dificuldade com condições de contorno • Velocidades no espaço paramétrico • Maior complexidade matemática • Malha deslocada: • Velocidades normais são mais intuitivas • Condições de contorno mais simples Velocidades paramétricas Malha colocalizada Malha deslocada

16. Advecção semi-Lagrangiana • Dados: • Tempo atual t, passo de integração h • Partícula numa posição amostral xAcom velocidade uA • Procedimento: • Partícula que seria carregada pelo fluido até o ponto amostral xA • Rastrear para trás (qualquer ordem, passo -h) até origem xB • Atribuir à amostra a velocidade interpolada em xB • Dissipação numérica • Interpolação cúbica saturada e outras técnicas • Efeito comparável à difusão/viscosidade

17. Advecção semi-Lagrangiana • Malhas curvilíneas • Integração pode ser feita no espaço global ou paramétrico • Maior eficiência no espaço em que a velocidade está representada • Maior precisão no espaço em que a velocidade exibe comportamento mais simples • Subdivisão em passos adaptativos • Ao atravessar fronteiras entre células ou percorrer distância unitária no espaço paramétrico • Considera variações nas velocidades e nos jacobianos • Atribuição final deve ser feita sempre com base no espaço do mundo • Evitar influência da variação de tamanho e rotações das células Trajetórias de uma partícula com velocidade u: A – linear no espaço do mundo (correta para u constante no mundo) B – linear no espaço paramétrico local (correta para u paramétrica constante) C – linear no espaço local com passo adaptativo (aproximada para u constante no mundo)

18. Difusão • Equação de Poisson implícita da difusão: [Stam99] • Malhas regulares (cada componente) • Solução de um sistema esparso simétrico com 7 diagonais em 3D (5 em 2D) • Iterações de Jacobi: Regra de atualização acessa 6 vizinhos em 3D (4 em 2D) • Exemplo de solver: Gradiente Conjugado (CG) com pré-condicionadorCholesky Incompleto (IC) • Malhas curvilíneas (cada componente) • Solução de um sistema esparso não-simétrico com 19 diagonais em 3D (9 em 2D) • Iterações de Jacobi: Regra de atualização acessa 18 vizinhos em 3D (8 em 2D) • Exemplo de solver: Gradiente Biconjugado Estabilizado (BiCGStab) com pré-condicionadorMultigridAlgébrico (AMG)

19. Forças externas • Integração explícita de primeira ordem é suficiente em qualquer caso • Acelerações arbitrárias violam fortemente condição de incompressibilidade • Aplicar projeção logo após esta etapa

20. Projeção • Visão geral • Após etapas anteriores, velocidades violam condição de incompressibilidade • Queremos projetar as velocidades em um campo livre de divergência • Utilizamos a decomposição de HelmHoltz-Hodge: • Procedimento • Calcular a divergência do campo de velocidades intermediárias • Calcular a pressão q (a menos de uma constante) resolvendo a equação de Poisson • Projetar a velocidade num campo livre de divergência: • Malhas regulares • Solução de um sistema esparso com a mesma estrutura do operador de difusão • Exemplo de solver: Gradiente Conjugado (CG) com pré-condicionador Cholesky Incompleto (IC)

21. Projeção em malhas curvilíneas • Divergência • Malha regular: • Malha curvilínea: • Espaço paramétrico: • As divergências das bases do espaço paramétrico podem ser pré-calculadas

22. Projeção em malhas curvilíneas • Pressão • Solução de um sistema esparso com a mesma estrutura do operador de difusão • Exemplo de solver: Gradiente Biconjugado Estabilizado (BiCGStab) com pré-condicionador Multigrid Algébrico (AMG) • Projeção • Malha regular: • Malha curvilínea: • Espaço paramétrico: • Os termos da matriz composta acima podem ser pré-calculados

23. Transporte de partículas e densidades de tinta • Transporte de partículas sem massa • Integração direta da posição das partículas • Segue as mesmas regras da advecção semi-lagrangeana • Posições guardadas diretamente no espaço paramétrico • Conversão do espaço do mundo para o paramétrico seria custosa • Advecção de densidades de tinta[Stam99] • Estrutura análoga à equação de Navier-Stokes • Operadores: advecção, difusão, extinção e fontes externas • Resolução análoga ao procedimento aplicado para Navier-Stokes

24. Condições de Contorno em malha regular • Camada extra de células junto às fronteiras e obstáculos • Cada célula de borda B está associada a um offset em (s,t,p) para uma célula interior R de referência • Relação entre pressões e velocidades de B e R determina diferentes comportamentos

25. Condições de Contorno em malha curvilínea • Aplicadas no espaço paramétrico normalizado • Componentes tangente e ortogonal sempre alinhadas com os eixos locais • Módulos independentes do tamanho local das células • Utilização de fatores de escala e offset: • Pressão: • Velocidade: • Malha deslocada e velocidades globais: • paredes lisas não suportadas devido àdificuldade em garantir que apenas acomponente normal seja igual a zero

26. Implementação – Dados CUDA • Propriedades que requerem leitura e escrita simultânea precisam de 2 áreas de memória (x2) • Cuda Arrays oferecem interpolação eficiente por hardware • Suporte a CUDA Surfaces depende do hardware do cliente • Dados com 3 componentes devem ser separados para manter o acesso coalescente

27. Implementação – Kernels CUDA • Não há necessidade de uso da memória compartilhada, pois todos os acessos são naturalmente coalescentes ou apresentam localidade espacial. • A ser testado: tornar os kernels limitados por computação em vez de por banda [Nguyen10] • Atualizações de células de borda: marcações “v” e “p” no diagrama • Gargalo: resolução de sistemas lineares esparsos (tarjas no diagrama) • Iterações de Jacobi • BiCGStab com pré-condicionador Multigrid Algébrico (implementados na biblioteca CUSP) • Visualização feita com o auxílio de shaders GLSL

28. Resultados – Validação Inicial • Resultados idênticos ao Stable Fluids em grades regulares 2D/3D

29. Resultados – Caminho curvo 2D

30. Resultados – Caminho curvo 3D

31. Resultados – Constrição 2D

32. Resultados – Constrição 3D

33. Resultados – Constrição 3D • Gráfico mostrando a velocidade média do fluido na direção x ao longo de seções transversais do grid com constrição. A área da constrição é 1/4 da área nas outras partes do percurso.

34. Resultados – Caminho suave 2D

35. Resultados – Caminho suave 3D

36. Resultados – Obstáculo 3D

37. Resultados • Informações gerais: • Intel Core i5 750 2.67GHz, 4GB RAM, placa gráfica GeForce GTX 550 Ti. • Gargalo: resolução das equações de Poisson • Convergência: Jacobi x BiCGStab • Iterações de Jacobi não conseguem atingir erros tão baixos quanto BiCGStab • BiCGStab não converge em alguns casos! • Melhora com uso da última pressão como chute inicial • Melhora significativamente com utilização de um ponto de pressão fixa

38. Resultados - Desempenho Tempos para a simulação completa (fluido, partículas e tinta), mas sem renderização 2 a 3 mais lento, mas melhor representação de domínios curvos

39. Conclusão • Método eficiente para a simulação euleriana de fluidos confinados em domínios tridimensionais curvilíneos estruturados • Matrizes Jacobianas relacionam o espaço do mundo com o paramétrico de malhas estruturadas • Simulação realizada com base numa simples grade uniforme • Desempenho • Técnica implementada em CUDA, explora bem o paralelismo das placas gráficas atuais. • 2 a 3 vezes mais lento que o Stable Fluids original com implementação similar. • Fronteiras curvilíneas • Dispensa informações topológicas • Não exige tratamento de casos especiais • Não requer maior refinamento próximo às fronteiras • Evita o surgimento de artefatos na simulação • Resolvedores de sistemas lineares esparsos • Iterações de Jacobi: mais simples de se implementar e paralelizar, estabilidade garantida • BiCGStab: convergência muito rápida, diminuindo o tempo de simulação e melhorando a qualidade do resultado, mas estabilidade não garantida • Trabalhos futuros • Abordagem integral para pressões • Técnicas de preservação da vorticidade • Tratamento de gravidade e superfície livre

40. Perguntas?

41. Referências [1] N Foster, D Metaxas, Modeling the Motion of a Hot, Turbulent Gas - Proceedings of ACM SIGGRAPH 1997 [2] J Stam, Stable Fluids – Proceedings of ACM SIGGRAPH 1999 [3] M J Harris, Fast Fluid Dynamics Simulation on the GPU – GPU Gems, 2004 [4] P Mullen et al., Energy-Preserving Integrators for Fluid Animation – Proceedings of ACM SIGGRAPH 2009 [5] E Wu, Y Liu, X Liu, An Improved Study of Real-Time Fluid Simulation on GPU – The Journal of Visualization and Computer Animation, 2004 [6] N Foster, R Fedkiw, Practical Animation of Liquids – Proceedings of ACM SIGGRAPH 2001 [7] D Enright, S Marschner, R Fedkiw, Animation and Rendering of Complex Water Surfaces – Proceedings of ACM SIGGRAPH 2002 [8] G Irving et al., Efficient Simulation of Large Bodies of Water by Coupling Two and Three Dimensional Techniques – Proceedings of ACM SIGGRAPH 2006 [9] J-M Hong, C-H Kim, Discontinuous Fluids – Proceedings of ACM SIGGRAPH 2005 [10] F Losasso et al., Multiple Interacting Liquids – Proceedings of ACM SIGGRAPH 2006 [11] A Robinson-Mosher et al., Two-way Coupling of Fluids to Rigid and Deformable Solids and Shells – Proceedings of ACM SIGGRAPH 2008 [12] E. Guendelman et al., Coupling Water and Smoke to Thin Deformable and Rigid Shells – Proceedings of ACM SIGGRAPH 2005 [13] C Batty, F Bertails, R Bridson, A Fast Variational Framework for Accurate Solid-Fluid Coupling - Proceedings of ACM SIGGRAPH 2007 [14] Nguyen A. et al., 3.5-D Blocking Optimization for Stencil Computations on Modern CPUs and GPUs –Proceedings of International Conference for High Performance Computing, Networking, Storage and Analysis, 2010 [15] AZEVEDO, V. C., Efficient smoke simulation on curvilinear grids. – Master’s thesis, UFRGS, 2012.

42. Anexos

43. Difusão em malhas regulares • Equação de Poisson implícita da difusão: [Stam99] • Malhas regulares:

44. Difusão em malhas curvilíneas