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
