wikimechanica.eesc.usp.br

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:

irrr = 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);

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.