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!

7 comments:

  1. excelente post!

    agora... digamos que eu quisesse alterar o laplacianFoam para resolver um problema diferente (tau d²T/dt² + dT/dt = alpha d²T/dx², conforme esse link); bastaria eu alterar o solve() no código fonte para

    solve (
    TAU * (fvm::d2dt2(T))
    + fvm::ddt(T)
    - fvm::laplacian(ALPHA, T)
    )

    , alterar também o 'createFields.H' definindo esses novos escalares (TAU, ALPHA) para defini-los depois no './constant/transportProperties', compilar e ser feliz, ou estou esquecendo de algo?

    abraço!

    ReplyDelete
  2. Obrigado pelo excelente post! Será muito útil para quem quiser aprofundar os seus conhecimentos do OpenFOAM.

    Um abraço,
    José Santos

    ReplyDelete
  3. @Anônimo: E que seja feliz!! A única coisa que eu faria é deixar o alpha fora do laplaciano, como:

    alpha*fvm::laplacian(T)

    Afinal de contas, é uma constante, não? Lembre-se de definir as constantes com as dimensões corretas. Depois é só correr pro abraço, pois tá prestes a fazer um gol!! Depois eu vou escrever sobre esse problema.

    @José Santos: Valeu pelo comentário, meu caro! Espero que seja de ajuda!

    Abraço!

    ReplyDelete
  4. hmmm.... olha só que interessante...
    eu pensei nisso também, mas na definição de laplacian() que vi tinham 2 argumentos, e tanto laplacian(1,T) quanto laplacian(1.,T) davam erro (1 != scalar); a maneira mais acéfala que encontrei foi colocar ALPHA dentro; :D

    ReplyDelete
  5. Anônimo,

    Para falar a verdade, a sua abordagem é válida sim e, numericamente, não vai fazer diferença no seu caso.

    Mas na programação da construção do fvm::laplacian é permitido usar as duas formas que discutimos. Veja no guia do programador do OpenFOAM.

    Abraço e boa sorte!

    ReplyDelete
  6. é, eu sei que é válida; mas, em contrapartida, usando laplacian(T) deve dar alguma diferença em termos de performance;

    aliás, falando nisso, e quanto a solvers 1D/2D? lembro de ter lido no programmer's guide que os operadores sempre assumiam a tridimensionalidade do domínio, mas gostaria de saber se alguém aqui poderia confirmar (ou negar) isso;

    ReplyDelete
  7. De fato, é uma questão de performance.

    O OpenFOAM realmente assume a tridimensionalidade do domínio. Contudo, vc pode ir eliminando uma (ou mais) dimensão (X, Y ou Z) definindo sua condição de contorno como "empty". O OpenFOAM não tem problema algum com geometrias 2D ou 1D. Leia o manual pois contém informações sobre isso.

    Além disso, veja o caso tutorial cavity, onde a dimensão Z é eliminada usando o contorno como empty. Dessa forma, a simulação é 2D.

    Um abraço.

    ReplyDelete