Nos testes realizados mais recentemente, observou-se que os solvers iterativos (GMRES, BiCGSTAB etc) não funcionam apropriadamente com elementos finitos triangulares quadráticos. No entanto, para elementos finitos lineares, são obtidos bons resultados. Os testes foram realizados no programa de escoamentos incompressíveis otimizado (matriz em banda com nós renumerados) para o escoamento sobre o aerofólio NACA 0012 com Re=1000. Os melhores resultados foram obtidos com o seguinte input: mpirun -np 4 ./f -pc_type bjacobi -ksp_type fgmres -ksp_gmres_restart 100 -ksp_max_it 100 -snes_monitor -snes_type newtonls -snes_linesearch_type l2 -snes_max_linear_solve_fail 3 -snes_max_it 3.

Observações:

Pré-condicionador: block Jacobi.

Solver: Tanto o GMRES quanto suas variantes (FGMRES, DGMRES etc) não pareceram apresentar muita diferença na convergência. O método BiCGSTAB também funcionou bem, inclusive sendo um pouco mais rápido do que o GMRES.

SNES: O problema foi resolvido dentro do paradigma SNES do PETSc, bastante recomendado pois, embora aumente levemente o tempo de processamento devido à integração separada da matriz jacobiana e do vetor RHS, proporciona uma convergência muito maior do que o Newton-Raphson puro (especialmente para o campo de pressão).

Problema: -Hélice 2D

ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT, 500);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc);
ierr = PCSetType(pc,PCASM);
ierr = KSPSetType(ksp,KSPGMRES); CHKERRQ(ierr);
ierr = KSPGMRESSetRestart(ksp, 500); CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);

Não converge a parte nao linear. O solver linear fica com erros relativamente pequenos.

Aumentando para 5000 iterações e 5 GMRES restart o não linear converge mas de forma diferente a cada nova simulação e o solver linear fica com erros grandes.

Mantendo 5000 iterações e aumentando para 50 GMRES o erro do solver linear diminui um pouco mas o não linear não converge bem, mas se ajusta com o passar do tempo.

Mantendo 5000 iterações e aumentando para 500 GMRES o erro do solver linear diminui mas a convergencia ainda é ruim, devendo se acertar com o processo dinâmico, eu espero.

Para o problema da hélice 2d resolvendo pelo método Arlequin o solver DGMRES sem precondicionamento ficou excelente. Converge em cerca de 30 a 90 iterações (um pouco mais nos primeiros passos de tempo), porém demanda um limite maior de iterações dentro do Newton-Raphson (utilizei em torno de 5).

ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-10,PETSC_DEFAULT,PETSC_DEFAULT, 5000);
CHKERRQ(ierr); ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPGetPC(ksp,&pc); ierr = PCSetType(pc,PCNONE);
ierr = KSPSetType(ksp,KSPDGMRES); CHKERRQ(ierr);
ierr = KSPGMRESSetRestart(ksp, 1000); CHKERRQ(ierr);
 ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);

O MUMPS (MUltifrontal Massively Parallel sparse direct Solver) é um solver direto, porém paralelo e muito eficiente. Para utilizá-lo em conjunto com o PETSc, a opção '–download-mumps' deve ser ativada no momento da instalação (ver PETSc). Sua utilização pode ser feita da seguinte forma (com todas as opções default).

    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
#if defined(PETSC_HAVE_MUMPS)
    ierr = KSPSetType(ksp,KSPPREONLY);
    ierr = KSPGetPC(ksp,&pc);
    ierr = PCSetType(pc, PCLU);
#endif
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
    ierr = KSPSetUp(ksp);
    ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);

No entanto, alguns erros podem acontecer durante a execução do MUMPs. Para monitorálos é importante ler o seu manual. De um modo geral, as tags de erro são impressas na variável INFO(1). Porém, essa informação só é acessada através do PETSc se a matriz tiver sido fatorada. Para fatorar uma matriz, adicionam-se os seguintes comandos:

ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
ierr = PCFactorSetUpMatSolverType(pc);
ierr = PCFactorGetMatrix(pc,&F);

Lembrando que uma nova matriz F deve ser declarada no preâmbulo da função.

Assim, para acessar a flag INFO(1) utiliza-se o seguinte comando:

PetscInt info1;
ierr = MatMumpsGetInfo(F,1,&info1);

info1);

Se INFO(1) = 0, o problema será realizado sem erros. Para os demais valores de INFO(1), consultar o manual de instruções do MUMPS.

Em um determinado problema obteve-se, por exemplo, a flag INFO(1) = -9. Segundo o manual de instruções do MUMPS, isso indica que o vetor interno S utilizados para os cálculos do MUMPS é muito pequeno. Para corrigir isso, deve-se aumentar o valor do parâmetro ICNTL(14), que é feito da seguinte forma:

#if defined(PETSC_HAVE_MUMPS)
        ierr = KSPSetType(ksp,KSPPREONLY);
        ierr = KSPGetPC(ksp,&pc);
        ierr = PCSetType(pc, PCLU);
 
        ierr = PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
        PCFactorSetUpMatSolverType(pc);
        PCFactorGetMatrix(pc,&F);
 
        PetscInt ival,icntl;
        icntl = 14; ival = 60;
        MatMumpsSetIcntl(F,icntl,ival);
#endif
 
        ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
        ierr = KSPSetUp(ksp);
 
#if defined(PETSC_HAVE_MUMPS)
        PetscInt  info1,info2,icntl14;
 
        MatMumpsGetInfo(F,1,&info1);
        MatMumpsGetInfo(F,2,&info2);
        MatMumpsGetIcntl(F,14,&icntl14);
        if(rank == 0) std::cout <<" INFO(1) = " <<info1
                                <<" " <<info2 <<" " <<icntl14 <<std::endl;
#endif

Note que a última parte do código é construída para monitorar os valores de INFO(1), INFO(2) e ICNTL(14). É fortemente recomendado que pelo menos INFO(1) e INFO(2) sejam monitorados a cada chamada do MUMPS para facilitar o diagnóstico de possíveis bugs ou erros de execução.