[petsc-users] 3.5 -> 3.6 Change in MatOrdering Behavior

[Paul T. Bauman]

Greetings,

We have a test that has started failing upon switching from 3.5.4 to 3.6.0
(actually went straight to 3.6.3 but checked this is repeatable with
3.6.0). I've attached the matrix generated with -mat_view binary and a
small PETSc program that runs in serial that reproduces the behavior by
loading the matrix and solving a linear system (RHS doesn't matter here).
For context, this matrix is the Jacobian of a Taylor-Hood approximation of
the axisymmetric incompressible Navier-Stokes equations for flow between
concentric cylinders (for which there is an exact solution). The matrix is
for a two element case, hopefully small enough for debugging.

Using the following command line options with the test program works with
PETSc 3.5.4 and gives a NAN residual with PETSc 3.6.0:

PETSC_OPTIONS="-pc_type asm -pc_asm_overlap 12 -sub_pc_type ilu
-sub_pc_factor_mat_ordering_type 1wd -sub_pc_factor_levels 4"

If I remove the mat ordering option, all is well again in PETSc 3.6.x:

PETSC_OPTIONS="-pc_type asm -pc_asm_overlap 12 -sub_pc_type ilu
-sub_pc_factor_levels 4"

Those options are nothing special. They were arrived at through trial/error
to get decent behavior for the solver on up to 4 processors to keep the
time to something reasonable for the test suite without getting really
fancy. Specifically, we'd noticed this mat ordering on some problems in the
test suite behaved noticeably better.

As always, thanks for your time.

Best,

Paul

test.mat
Description: Binary data

test.mat.info
Description: Binary data

#include <petscksp.h>

int main(int argc,char **args)
{
  Vec            u,b;
  Mat            A;
  KSP            ksp;
  PetscViewer    viewer;
  PetscInt       n_rows, n_cols;
  PetscErrorCode ierr;

  char filename[] = "./test.mat";

  PetscInitialize(&argc,&args,(char*)0,"");

  ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);

  ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
  ierr = MatSetFromOptions(A);CHKERRQ(ierr);
  ierr = MatLoad(A,viewer);CHKERRQ(ierr);

  ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

  ierr = MatGetSize(A,&n_rows,&n_cols);CHKERRQ(ierr);

  ierr = VecCreate(PETSC_COMM_WORLD,&b);CHKERRQ(ierr);
  ierr = VecSetSizes(b,n_rows,PETSC_DECIDE);CHKERRQ(ierr);
  ierr = VecSetFromOptions(b);CHKERRQ(ierr);

  ierr = VecDuplicate(b,&u);CHKERRQ(ierr);

  ierr = VecSet(b,1.0);CHKERRQ(ierr);

  ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
  ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
  ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
  ierr = KSPSolve(ksp,b,u);CHKERRQ(ierr);

  ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
  ierr = MatDestroy(&A);CHKERRQ(ierr);

  ierr = PetscFinalize();

  return 0;
}