[firedrake] Extracting components of FieldSplit preconditioner

[Eike Mueller]

Dear firedrakers,

I’m trying to extract the components of the field split preconditioner to time them individually, but it fails. I suspect that I’mnot setting up some things, but it’s not obvious what to do. For example, if I do 

        ksp = up_solver.snes.getKSP()
        ksp.setUp()
        ksp_hdiv, ksp_schur = ksp.getPC().getFieldSplitSubKSP()

        # HDiv space
        ksp_hdiv.setUp()
        op_hdiv, op_pc_hdiv = ksp_hdiv.getOperators()
        pc_hdiv = ksp_hdiv.getPC()
        print pc_hdiv.getType()
        ndof = op_hdiv.getSizes()[0][0]
        x = PETSc.Vec()
        x.create()
        x.setSizes((ndof, None))
        x.setFromOptions()
        x.setArray(np.random.rand(ndof))
        y = x.duplicate()
        x.setArray(np.zeros(ndof))
        pc_hdiv.apply(x,y)
        with timed_region('pc_hdiv'):
            pc_hdiv.apply(x,y)

I get:

bjacobi
Traceback (most recent call last):
  File "driver.py", line 557, in <module>
    main(parameter_filename)
  File "driver.py", line 499, in main
    pc_hdiv.apply(x,y)
  File "PC.pyx", line 187, in petsc4py.PETSc.PC.apply (src/petsc4py.PETSc.c:132142)
petsc4py.PETSc.Error: error code 71
[0] PCApply() line 449 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/interface/precon.c
[0] PCApply_BJacobi_Singleblock() line 670 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/impls/bjacobi/bjacobi.c
[0] KSPSolve() line 542 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/ksp/interface/itfunc.c
[0] KSPSetUp() line 330 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/ksp/interface/itfunc.c
[0] PCSetUp() line 918 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/interface/precon.c
[0] PCSetUp_ILU() line 232 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/ksp/pc/impls/factor/ilu/ilu.c
[0] MatLUFactorNumeric() line 2945 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/mat/interface/matrix.c
[0] MatLUFactorNumeric_SeqAIJ() line 566 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/src/mat/impls/aij/seq/aijfact.c
[0] MatPivotCheck() line 664 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/include/petsc-private/matimpl.h
[0] MatPivotCheck_none() line 645 in /Users/eikemueller/PostDocBath/EllipticSolvers/petsc/include/petsc-private/matimpl.h
[0] Zero pivot in LU factorization: http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot
[0] Zero pivot row 0 value 0 tolerance 2.22045e-14

so it looks to me like the PC is not set up correctly.

I use

                 'fieldsplit_0_ksp_type': 'preonly',
                 'fieldsplit_0_pc_type': 'bjacobi',
                 'fieldsplit_0_sub_pc_type': 'ilu’,

So I suspect the ILU has not been set up correctly

Thanks,

Eike

The solver parameters are:

 sparams={'pc_type': 'fieldsplit',
                 'pc_fieldsplit_type': 'schur',
                 'ksp_type': 'gmres',
                 'ksp_max_it': 100,
                 'ksp_rtol':self._tolerance,
                 'pc_fieldsplit_schur_fact_type': 'FULL',
                 'pc_fieldsplit_schur_precondition': 'selfp',
                 'fieldsplit_0_ksp_type': 'preonly',
                 'fieldsplit_0_pc_type': 'bjacobi',
                 'fieldsplit_0_sub_pc_type': 'ilu',
                 'fieldsplit_1_ksp_type': 'preonly',
                 'fieldsplit_1_pc_type': 'gamg',
                 'fieldsplit_1_mg_levels_ksp_type': 'chebyshev',
                 'fieldsplit_1_mg_levels_ksp_chebyshev_estimate_eigenvalues': True,
                 'fieldsplit_1_mg_levels_ksp_chebyshev_estimate_eigenvalues_random': True,
                 'fieldsplit_1_mg_levels_ksp_max_it': 1,
                 'fieldsplit_1_mg_levels_pc_type': 'bjacobi',
                 'fieldsplit_1_mg_levels_sub_pc_type': 'ilu',
                 'ksp_monitor': False}