February 10, 2009

Análise de código no OpenFOAM - laplacianFoam

Para explicar a estrutura e o algoritmo de solução de um código no OpenFOAM, vou usar como base o solver laplacianFoam. Para acompanhar melhor esse post, é interessante que o leitor tenha alguns conhecimentos básicos da sintaxe de C++. Porém, vou apresentar alguns detalhes referentes aos comandos e funções que são membros das classes e templates, facilitando a leitura do código para os leitores sem experiência em linguagens orientadas a objetos.

O solver laplacianFoam é usado para resolver o problema da difusão pura de um campo escalar T, sem considerar nenhum termo fonte. Esta equação está colocada abaixo, sendo D o coeficiente de difusão.

Os arquivos referentes aos solvers do OpenFOAM ficam no diretório OpenFOAM-version/applications/solvers., onde version se refere a versão do OpenFOAM. O código do laplacianFoam fica no diretório basic e está colocado abaixo.
00001 /*---------------------------------------------------------------------------*\
00002 ========= |
00003 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
00004 \\ / O peration |
00005 \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
00006 \\/ M anipulation |
00007 -------------------------------------------------------------------------------
00008 License
00009 This file is part of OpenFOAM.
00010
00011 OpenFOAM is free software; you can redistribute it and/or modify it
00012 under the terms of the GNU General Public License as published by the
00013 Free Software Foundation; either version 2 of the License, or (at your
00014 option) any later version.
00015
00016 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
00019 for more details.
00020
00021 You should have received a copy of the GNU General Public License
00022 along with OpenFOAM; if not, write to the Free Software Foundation,
00023 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
00024
00025 Application
00026 laplacianFoam
00027
00028 Description
00029 Solves a simple Laplace equation, e.g. for thermal diffusion in a solid.
00030
00031 \*---------------------------------------------------------------------------*/
00032
00033 #include "fvCFD.H"
00034
00035
00036 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00037
00038 int main(int argc, char *argv[])
00039 {
00040
00041 # include "setRootCase.H"
00042
00043 # include "createTime.H"
00044 # include "createMesh.H"
00045 # include "createFields.H"
00046
00047 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
00048
00049 Info<< "\nCalculating temperature distribution\n" << endl;
00050
00051 for (runTime++; !runTime.end(); runTime++)
00052 {
00053 Info<< "Time = " << runTime.timeName() << nl << endl;
00055 # include "readSIMPLEControls.H"
00056
00057 for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
00058 {
00059 solve
00060 (
00061 fvm::ddt(T) - fvm::laplacian(DT, T)
00062 );
00063 }
00064
00065 # include "write.H"
00066
00067 Info<< "ExecutionTime = " << runTime.elapsedCpuTime() <<" s"
00068 << " ClockTime = " << runTime.elapsedClockTime() <<" s"
00069 << nl << endl;
00070 }
00071
00072 Info<< "End\n" << endl;
00074 return(0);
00075 }
00076
00077
00078 // ************************************************************************* //
Para acompanhar a leitura, sugiro que você abra diretamente o arquivo laplacianFoam.C em um editor de textos de sua preferência. De forma mais direta, você também acessar o código do laplacianFoam abrindo este link em uma outra aba do seu navegador.

A primeira linha a ser executada no código (linha 33) declara a biblioteca fvCFD.H, que fornece ao solver acesso a todas as classes e propriedades descritas em meu post anterior. Deve-se frisar que é vital declarar esta biblioteca. Na linha 38, a função main engloba todo o código fonte principal e possui dois argumentos de entrada: o inteiro argc e a cadeia de caracteres argv. Estes parâmetros contém informações sobre a simulação, como o diretório e o nome do caso a ser simulado. Os argumentos são lidos pelo programa diretamente na linha de comando para execução do solver.

A biblioteca setRootCase.H é usada para testar a validade dos argumentos argc e argv da simulação. O conteúdo desta biblioteca está colocado abaixo.
00001 Foam::argList args(argc, argv);
00002
00003 if (!args.checkRootCase())
00004 {
00005 Foam::FatalError.exit();
00006 }
O primeiro comando do código acima declara a variável args construída com os argumentos argc e argv, a partir da classe argList. Em seguida, o comando checkRootCase() verifica a validade e a existência do diretório e do nome do caso simulado. Caso o retorno de checkRootCase() seja False, a execução do solver é interrompida pelo comando padrão de erro do OpenFOAM FatalError.

As duas próximas bibliotecas declaradas, createTime.H e createMesh.H, são responsáveis pela criação de bancos de dados para armazenar dados sobre o caso simulado e a estrutura da malha utilizada. O código colocado abaixo refere-se à biblioteca createTime.H.
00001     Foam::Info<< "Create time\n" << Foam::endl;
00002
00003 Foam::Time runTime
00004 (
00005 Foam::Time::controlDictName,
00006 args.rootPath(),
00007 args.caseName()
00008 );
Sendo construída a partir de informações sobre o nome e o diretório do caso simulado (provindas de args), a variável runTime do template Time é declarada no código acima. Desta forma, runTime obtém a localização do arquivo de configuração controlDict do caso e, utilizando as informações contidas neste último, monta um banco de dados para controle da simulação. Por exemplo, pode-se citar alguns dos dados contidos em runTime: (i) instante inicial e final; (ii) controle do passo de tempo (fixo, adaptativo, etc.); (iii) diretórios que contém os arquivos com os campos iniciais das propriedades transportadas; (iv) controle de escrita em arquivo (formato de saída, compressão de dados, etc.); entre outros.

A biblioteca createMesh.H usa o template fvMesh para declarar a variável mesh, construída a partir de outro template chamado IOobject, como mostra o código a seguir.
00001     Foam::Info<< "Create mesh for time = "
00002 << runTime.timeName()<< Foam::nl << Foam::endl;

00003
00004 Foam::fvMesh mesh
00005 (
00006 Foam::IOobject
00007 (
00008 Foam::fvMesh::defaultRegion,
00009 runTime.timeName(),
00010 runTime,
00011 Foam::IOobject::MUST_READ
00012 )
00013 );
O template IOobject define os atributos de um objeto de modo a fornecer meios para entrada e/ou saída de dados (usualmente em arquivo). Com as informações de runTime, o template fvMesh é capaz de localizar os arquivos cells, faces, points e boundary para construção da malha. Note que o último parâmetro na construção de IOobject refere-se à regras de leitura e escrita de arquivos. Neste caso, a leitura dos dados deve ser obrigatória (MUST_READ).

Todas as bibliotecas supracitadas são gerais e podem ser usadas em qualquer código do OpenFOAM (salvo pequenas modificações já comentadas). Contudo, a biblioteca createFields.H é usada para criação e leitura dos campos iniciais das incógnitas do problema e leitura de propriedades físicas aplicadas a cada caso (propriedades de transporte, termodinâmicas, termofísicas, etc.). Desta forma, este header é específico para cada solver e deve ser desenvolvido com cuidado pelo programador, pois todas as incógnitas e todas as propriedades físicas do problema devem ser definidas neste arquivo. Portanto, deve-se ter domínio do modelo fluidodinâmico e qual a melhor forma de armazenar suas variáveis de modo a otimizar o código. A biblioteca createFields.H específica para o laplacianFoam está colocada abaixo.
00001     Info<< "Reading field T\n" << endl;
00002
00003 volScalarField T
00004 (
00005 IOobject
00006 (
00007 "T",
00008 runTime.timeName(),
00009 mesh,
00010 IOobject::MUST_READ,
00011 IOobject::AUTO_WRITE
00012 ),
00013 mesh
00014 );
00015
00016
00017 Info<< "Reading transportProperties\n" << endl;
00018
00019 IOdictionary transportProperties
00020 (
00021 IOobject
00022 (
00023 "transportProperties",
00024 runTime.constant(),
00025 mesh,
00026 IOobject::MUST_READ,
00027 IOobject::NO_WRITE
00028 )
00029 );
00030
00031
00032 Info<< "Reading diffusivity DT\n" << name="l00033">00033
00034 dimensionedScalar DT
00035 (
00036 transportProperties.lookup("DT")
00037 );
Para resolver o problema, a temperatura deve ser definida em um campo geométrico (template geometricField). Fisicamente, a temperatura é uma variável escalar e, portanto, deve-se criar um campo geométrico de um escalar. O template volScalarField constrói um campo escalar da variável T, definida no centro das células da malha mesh. A criação deste campo depende dos templates IOobject e fvMesh como argumentos de entrada. O primeiro constrói o objeto, definindo o nome da variável e do arquivo que contém os valores do campo inicial de T (O nome do arquivo que contém o campo inicial da variável é idêntico ao nome da própria variável), o registro das informações contidas em runTime (por exemplo, o diretório e a periodicidade da saída de resultados em arquivo) e as opções de entrada e saída de dados (MUST_READ indica que o campo inicial de T deve ser necessariamente lido e AUTO_WRITE configura a saída automática de T em arquivo ao longo do tempo). O segundo parâmetro para a criação de T é o template fvMesh definido como mesh, que indica onde o campo escalar T será alocado e inserido.

A equação da difusão possui apenas uma propriedade física de transporte, a condutividade térmica DT (representada pela letra D na equação da difusão). A leitura das propriedades de transporte é realizada através do template IOdictionary, construída a partir do template IOobject como argumento. O template IOdictionary, por sua vez, é derivada de outros dois templates, dictionary e IOobject, proporcionando funcionalidade na entrada e saída de dados automática a partir de um banco de dados. O template dictionary define uma lista de palavras chave, onde cada uma destas é associada a um número arbitrário de valores. A construção de transportProperties declara o nome do arquivo que contém as propriedades de transporte e seu diretório (runTime.constant()) e as regras de leitura (MUST_READ) e saída (NO_WRITE) de dados em arquivo. Por fim, a criação da variável escalar dimensional DT é construída pelo comando transportProperties.lookup("DT"), que procura no arquivo transportProperties a palavra-chave "DT" e associa um valor dimensional a essa variável.

As etapas descritas acima apenas inicializam os dados para runTime, a malha mesh, o campo inicial de T e a propriedade de transporte DT do código principal. O próximo passo é programar o algoritmo de solução do problema específico.

O laço for, iniciado na linha 51 do código principal, tem o intuito de repetir as instruções no interior do laço para cada passo de tempo (incrementado por runTime++). O laço é realizado até que a condição de runTime.end() seja satisfeita (retorne true).

O primeiro comando no interior do laço é a declaração do header readSIMPLEControls.H para leitura dos parâmetros do método SIMPLE de acoplamento pressão-velocidade e das condições de ortogonalidade da malha. Estes parâmetros são lidos no arquivo fvSolution do caso analisado. Apesar de não ser necessário aplicar o método SIMPLE para resolver a equação da difusão pura, esta biblioteca lê o número de iterações para correção dos fluxos devido à não-ortogonalidade da malha. Isso é necessário pois o OpenFOAM divide o cálculo do fluxo através das faces em duas parcelas chamadas contribuição ortogonal e a correção não ortogonal. Esta correção é realizada um laço para ajustar o fluxo das propriedades nas faces dos volumes, semelhante ao esquema defferred correction. Maiores detalhes podem ser encontrados na tese do Prof. Jasak.

A discretização por volumes finitos é realizada pelo template fvm, armazenando as equações discretizadas em sua forma matricial com o template fvMatrix e montando um sistema de equações lineares resolvido pelo comando solve. Este último comando retorna na tela para o usuário dados estatísticos da solução, como a convergência do sistema, número de iterações, etc. A definição das formulações de discretização usadas na simulação estão colocadas no arquivo fvShemes. A solução do sistema algébrico, cujo método de solução está selecionado no arquivo fvSolution, retorna a temperatura em cada célula da malha. A biblioteca write.H, colocada abaixo, escreve os arquivos de resultados quando o comando de classe runTime.output() for válido.
00001     if (runTime.outputTime())
00002 {
00003 volVectorField gradT = fvc::grad(T);
00004
00005 volScalarField gradTx
00006 (
00007 IOobject
00008 (
00009 "gradTx",
00010 runTime.timeName(),
00011 mesh,
00012 IOobject::NO_READ,
00013 IOobject::AUTO_WRITE
00014 ),
00015 gradT.component(vector::X)
00016 );
00017
00018 volScalarField gradTy
00019 (
00020 IOobject
00021 (
00022 "gradTy",
00023 runTime.timeName(),
00024 mesh,
00025 IOobject::NO_READ,
00026 IOobject::AUTO_WRITE
00027 ),
00028 gradT.component(vector::Y)
00029 );
00030
00031 volScalarField gradTz
00032 (
00033 IOobject
00034 (
00035 "gradTz",
00036 runTime.timeName(),
00037 mesh,
00038 IOobject::NO_READ,
00039 IOobject::AUTO_WRITE
00040 ),
00041 gradT.component(vector::Z)
00042 );
00043
00044
00045 runTime.write();
00046 }
Com o intuito de também escrever em arquivo os componentes do gradiente da temperatura, calcula-se gradT pela classe fvc que realiza operações tensoriais explícitas com os dados da malha. Note que, como a temperatura é uma variável escalar, seu gradiente será uma entidade vetorial, e portanto o tipo da variável gradT deve ser um volVectorField. As linhas 5, 18 e 31 do código acima decompõem gradT nas direções X, Y e Z do domínio e definem esta variável como AUTO_WRITE para que a saída em arquivo seja automática. Cada decomposição é, portanto, um volScalarField e será formado pela componente do vetor de gradT da respectiva direção. O comando runTime.write() escreve em arquivos os valores de temperatura nos centros das células, seus gradientes e os componentes do gradiente.

O algoritmo repete-se até o fim do laço no tempo e, por consequência, o final da simulação. Na minha opinião, o solver laplacianFoam possui o código mais simples de todos. A complexidade dos códigos aumenta junto com o detalhamento dos fenômenos físicos considerados, sendo necessário o desenvolvimento de algoritmos mais elaborados, o uso de mais templates, classes e comandos.

Espero que este post seja de ajuda aos que estão começando a estudar e desbravar a programação no OpenFOAM. Se quiserem discutir mais sobre o assunto ou mesmo trocar idéias, deixe um comentário aqui.

Um abraço!