[petsc-users] ILU and Block Gauss Seidel smoothing

[Barry Smith]

On May 25, 2011, at 2:45 AM, Gianluca Meneghello wrote:

> Dear all,
> 
> I am writing a multigrid solver for the Navier Stokes Equations, using
> PETSc. I am not using the dmmg structure because I need grid
> refinement and, as far as I can see, that is not supported.
> 
> I have two questions concerning two different possibilities of smoothing:
> 
> 1) ILU relaxation: in Multigrid (U. Trottenberg,Cornelis W.
> Oosterlee,Anton Schüller), it is suggested that alternating ILU
> represents a robust smoother. Given a routine for ILU decomposition,
> my understanding is that alternation depends on the ordering of the
> matrix. My code builds a linearized Navier Stokes operator where first
> the u momentum equation, then the v momentum and then the p equation
> are stored. Each equation is built in EN order — first comes the
> bottom horizontal line, from West to East, then I move line by line
> from South to North. Otherwise stated, the position of each unknown in
> the matrix is given by P = cf*nx*ny + j*nx + i, where cf is the
> variable number (0 = u, 1 = v, 2 = p) and i,j,nx,ny are as usual.
> 
> For "historical" reasons I'm not using the DM structure even if the
> grids are logically rectangular, but I can change that.
> 
> My question: is there in PETSc a way of doing alternating ILU which
> does not require rebuilding the matrix with different ordering each
> time?

    There is a way for you to set the ordering of the ILU factorization without building the matrix in that ordering. It requires that you write a function with the calling sequence

PetscErrorCode  YourOrdering(Mat mat,const MatOrderingType type,IS *isrow, IS *iscol) 
{
    Your function fills up isrow with the order of the rows you want the ILU to visit, for example if you want the first row visited to be ten then you would create an isrow whose first entry is 10. You can/and probably should use the same ordering for rows and columns. Your function will igore the type name that is passed in.
}

   In your main program after PetscInitialize() call MatOrderingRegister("yourname",0,yourname,YourOrdering); 

now when you run the program you can use the option -pc_factor_mat_ordering_type yourname  and it will use the function you provided to generate the ordering that is used for ILU.

   Good luck and sorry for not responding to this message sooner,

   Barry


> 
> Somehow related questions: what is the best ordering for the matrix
> and which one is the one used by the DM structure?
> 
> 
> 2) Decoupled block-line Gauss Seidel relaxation: this is another
> smoother I'm considering, which explains the matrix ordering above (I
> use block horizontal lines smoothing, so that each line is stored
> consecutively in the matrix). At the moment the full linearized
> operator is first built and then parts of the matrix are extracted.
> In reality though, I never need the full linearized operator to be
> built except on the coarsest grid. Rather, I need it to be constructed
> each block-line at a time. I've seen that for parallel applications
> each part of the matrix is built on its corresponding processor. Is
> there a way to do it sequentially and to have control on which part is
> built?
> 
> 
> Any other suggestion on both subjects is of course welcome. Thanks in advance
> 
> Gianluca
> 
> 
> 
> -- 
> "[Je pense que] l'homme est un monde qui vaut des fois les mondes et
> que les plus ardentes ambitions sont celles qui ont eu l'orgueil de
> l'Anonymat" -- Non omnibus, sed mihi et tibi
> Amedeo Modigliani