[petsc-users] ILU and Block Gauss Seidel smoothing

[Gianluca Meneghello]

Barry,

thanks for your answer and sorry for my late reply.

While your suggestion is in the right direction for what I need, I
think I would need something a little different.

I would like to use different ordering in the same run, corresponding
to the different ordering of the unknown in a rectangular grid — that
is in NW, NE, SW, SE directions.
My goal is to perform some relaxation sweeps for each ordering at each
level of a multigrid process, probably using ksptype PREONLY and
pctype ILU.
Is that possible? Is there the equivalent of -pc_sor_its with ILU
(-pc_ilu_its maybe)?

I also have the problem that in order to build the ordering I would
need to have access to a structure containing some grid informations,
and it seems I cannot pass that structure to the YourOrdering function
you suggested.

I guess a solution for me could be to build the IS from an external
function (not used by petsc) and then attach them directly to the mat
structure. I also guess that the one to use are the ones at line 36 of

http://www.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-dev/src/mat/impls/aij/seq/aij.h.html

Not being an expert in the underlying structure of PETSc, I wonder
whether this is a good idea or not. Is there a more standard
alternative?

Thanks

Gianluca

On 3 June 2011 17:20, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> 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
>
>



-- 
"[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