[petsc-users] ILU and Block Gauss Seidel smoothing

[Barry Smith]

  richardson, you may also need to set that KSP with KSPSetInitialGuessNonzero() 

   Barry



On Jul 18, 2011, at 1:36 PM, Gianluca Meneghello wrote:

> Barry, 
> I'm actually using ilu as a smoother inside a multigrid cycle. Does in this case make sense to use ksppreonly? Or should I use Richardson?
> 
> Thanks
> 
> 
> 
> Il giorno 18/lug/2011, alle ore 20:14, Barry Smith <bsmith at mcs.anl.gov> ha scritto:
> 
>> 
>> On Jul 18, 2011, at 9:43 AM, Gianluca Meneghello wrote:
>> 
>>> Barry,
>>> 
>>> thanks again for your input. That would be really good.
>>> 
>>> When executing my program with the options you suggested I obtain
>>> various errors (unknown pc etc) so I tried to figure out how to avoid
>>> them. Is it something like this what you intended?
>>> 
>>> ./main
>>> -mg_pc_type composite
>>> -mg_pc_composite_pcs ilu,ilu
>>> -mg_ksp_type preonly
>> 
>>   I don't think you want this one. It says not to use a Krylov method on the outside; the preconditioner will just apply the two ILU triangular and that is it.  In other words it won't even solve the linear system. You likely want -mg_ksp_type gmres
>> 
>>  Barry
>> 
>>> -mg_sub_0_pc_factor_mat_ordering_type FESOrder
>>> -mg_sub_1_pc_factor_mat_ordering_type FWSOrder
>>> -mg_pc_composite_type multiplicative
>>> -mg_ksp_view
>>> 
>>> provided that
>>> 
>>> KSPSetOptionsPrefix(ksp,"mg_");
>>> 
>>> has been applied to the KSP I want to modify.
>>> 
>>> Thanks again
>>> 
>>> Gianluca
>>> 
>>> 
>>> On 13 July 2011 15:33, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>> 
>>>> On Jul 13, 2011, at 4:02 AM, Gianluca Meneghello wrote:
>>>> 
>>>>> Thanks Barry.
>>>>> 
>>>>> What I'm doing at the moment in order to pass the ordering to the PC is
>>>>> 
>>>>> PCFactorSetMatOrderingType(PC,"FESOrder");
>>>>> 
>>>>> where the FESOrder(...) function computes the ordering and has been
>>>>> registered in the main function with
>>>>> 
>>>>> MatOrderingRegister("FESOrder",0,"FESOrder",FESOrder);
>>>>> 
>>>>> The same is done for al the different orderings I'm using inside the
>>>>> PCSHELL suggested by Jed.
>>>>> 
>>>>> My understanding is that the FESOrder(...) function is evaluated every
>>>>> time the PCApply(...) function is called.
>>>> 
>>>> Definitely not. Looking at PCSetUp_ILU()
>>>> 
>>>> } else {
>>>>  if (!pc->setupcalled) {
>>>>    /* first time in so compute reordering and symbolic factorization */
>>>>    ierr = MatGetOrdering(pc->pmat,((PC_Factor*)ilu)->ordering,&ilu->row,&ilu->col);CHKERRQ(ierr);
>>>> 
>>>> it will actually call the ordering routine a TOTAL of ONE TIME if you use SAME_NONZERO_PATTERN to KSPSetOperators() or PCSetOperators(). (it is smart enough to reuse the ordering).
>>>> 
>>>> BTW: I don't think you even need to use a PCSHELL preconditioner. You can get the same effect with a PCCOMPOSITE.  Just create a single KSP in the usual way and use the options
>>>> 
>>>> -pc_type composite -pc_composite_type ilu,ilu -sub_0_ksp_type preonly -sub_1_ksp_type preonly -sub_0_pc_factor_mat_ordering_type FESOrder -sub_1_pc_factor_mat_ordering_type anotherorder somethingelse
>>>> 
>>>> 
>>>> Barry
>>>> 
>>>> 
>>>> 
>>>> 
>>>> 
>>>>> 
>>>>> If there was a way to store the IS containing the ordering — or any
>>>>> other array — into the shell context (this I know how to do that) and
>>>>> the pass the IS directly to the preconditioner instead of using
>>>>> PCFactorSetMatOrderingType would be the best solution for me.
>>>>> 
>>>>> Thanks
>>>>> 
>>>>> Gianluca
>>>>> 
>>>>> 
>>>>> On 12 July 2011 20:15, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>>>> 
>>>>>> On Jul 12, 2011, at 3:47 AM, Gianluca Meneghello wrote:
>>>>>> 
>>>>>>> Jed, thanks for your answer.
>>>>>>> 
>>>>>>> This is exactly what I was looking for (it took me some time to test it...).
>>>>>>> 
>>>>>>> Let me ask you a couple of additional questions:
>>>>>>> 
>>>>>>> 1) in this configuration, is there an alternative way to pass the
>>>>>>> ordering to the preconditioner? Becuse my ordering depends on the grid
>>>>>>> an not on the solution, it is constant during all the computations and
>>>>>>> I have no advantage on calling the function that creates multiple
>>>>>>> times (the pc is built and destroyed several times).
>>>>>> 
>>>>>> Not sure exactly what you want here. You can put the ordering information in the shell context and just keep it there forever and don't need to to recompute.
>>>>>> 
>>>>>>  PCShellSet/GetContext()?
>>>>>>> 
>>>>>>> 2) what's the best way to make the shell pc selectable from the option
>>>>>>> database, so that it can be replaced with more standard ones?
>>>>>>> 
>>>>>> 
>>>>>>      /* here the user can provide -pc_type shell or ilu or whatever */
>>>>>>      KSPSetFromOptions(ksp);
>>>>>>      KSPGetPC(ksp,&pc);
>>>>>>      /* the following lines will only do something if the user selected the PCType of shell otherwise they are ignored */
>>>>>> 
>>>>>>      PCShellSetContext(pc,....)
>>>>>>      PCShellSetSetUp(pc,....)
>>>>>>      PCShellSetApply(pc,.....)
>>>>>> 
>>>>>> 
>>>>>> Barry
>>>>>> 
>>>>>> 
>>>>>>> Thanks
>>>>>>> 
>>>>>>> Gianluca
>>>>>>> 
>>>>>>> On 1 July 2011 14:43, Jed Brown <jed at 59a2.org> wrote:
>>>>>>>> On Fri, Jul 1, 2011 at 07:11, Gianluca Meneghello <gianmail at gmail.com>
>>>>>>>> wrote:
>>>>>>>>> 
>>>>>>>>> I've tried to implement Jed's solution as it was looking faster to
>>>>>>>>> try. This is what I came up with for two different orders, FENOrder
>>>>>>>>> and FWNOrder respectively
>>>>>>>>> 
>>>>>>>>> 283       KSP ksp;
>>>>>>>>> 284       PC pc;
>>>>>>>>> 285
>>>>>>>>> 286       ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);  CHKERRQ(ierr);
>>>>>>>>> 287       ierr = KSPSetOptionsPrefix(ksp,"mg_");    CHKERRQ(ierr);
>>>>>>>>> 288       ierr = KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN);
>>>>>>>>> CHKERRQ(ierr);
>>>>>>>>> 289       ierr = KSPSetFromOptions(ksp);    CHKERRQ(ierr);
>>>>>>>>> 290       ierr = KSPGetPC(ksp,&pc);         CHKERRQ(ierr);
>>>>>>>>> 291       ierr = PCFactorSetFill(pc,1);     CHKERRQ(ierr);
>>>>>>>>> 292       ierr = PCFactorSetMatOrderingType(pc,"FENOrder"); CHKERRQ(ierr);
>>>>>>>>> 293       ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);
>>>>>>>>> 294       ierr = KSPDestroy(&ksp);  CHKERRQ(ierr);
>>>>>>>>> 295
>>>>>>>>> 296       ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);  CHKERRQ(ierr);
>>>>>>>>> 297       ierr = KSPSetOptionsPrefix(ksp,"mg_");    CHKERRQ(ierr);
>>>>>>>>> 298       ierr = KSPSetOperators(ksp,J,J,SAME_NONZERO_PATTERN);
>>>>>>>>> CHKERRQ(ierr);
>>>>>>>>> 299       ierr = KSPSetFromOptions(ksp);    CHKERRQ(ierr);
>>>>>>>>> 300       ierr = KSPGetPC(ksp,&pc);         CHKERRQ(ierr);
>>>>>>>>> 301       ierr = PCFactorSetFill(pc,1);     CHKERRQ(ierr);
>>>>>>>>> 302       ierr = PCFactorSetMatOrderingType(pc,"FWNOrder"); CHKERRQ(ierr);
>>>>>>>>> 303       ierr = KSPSolve(ksp,b,x); CHKERRQ(ierr);
>>>>>>>>> 304       ierr = KSPDestroy(&ksp);  CHKERRQ(ierr);
>>>>>>>> 
>>>>>>>> It looks like you are putting all of this in PCApply_YourShell(). I suggest
>>>>>>>> making
>>>>>>>> struct YourShellContext {
>>>>>>>> PC fenpc,fwnpc;
>>>>>>>> Vec work;
>>>>>>>> PetscBool issetup;
>>>>>>>> };
>>>>>>>> Then pass a pointer to this thing in when you create the PCShell. Then
>>>>>>>> something like
>>>>>>>> PCSetUp_YourShell(PC pc) {
>>>>>>>> struct YourShellContext *ctx;
>>>>>>>> Mat A,B;
>>>>>>>> MatStructure mstr;
>>>>>>>> const char *prefix;
>>>>>>>> char iprefix[256];
>>>>>>>> PCShellGetContext(pc,&ctx);
>>>>>>>> PCGetOptionsPrefix(pc,&prefix);
>>>>>>>> if (!ctx->issetup) {
>>>>>>>> ierr = MatGetVecs(A,&ctx->work,PETSC_NULL);
>>>>>>>> PCCreate(comm,&ctx->fenpc);
>>>>>>>> PCSetType(ctx->fenpc,PCILU);
>>>>>>>> PCFactorSetFill(ctx->fenpc,1);
>>>>>>>> PCFactorSetMatOrderingType(ctx->fenpc,"FENOrder");
>>>>>>>> snprintf(iprefix,sizeof iprefix,"%sfen",prefix);
>>>>>>>> PCSetOptionsPrefix(ctx->fenpc,iprefix);
>>>>>>>> PCSetFromOptions(ctx->fenpc);
>>>>>>>> /* Same as above for ctx->fwnpc */
>>>>>>>> }
>>>>>>>> PCGetOperators(pc,&A,&B,&mstr);
>>>>>>>> PCSetOperators(ctx->fenpc,A,B,mstr);
>>>>>>>> PCSetOperators(ctx->fwnpc,A,B,mstr);
>>>>>>>> PCSetUp(ctx->fenpc);
>>>>>>>> PCSetUp(ctx->fwnpc);
>>>>>>>> ctx->issetup = PETSC_TRUE;
>>>>>>>> }
>>>>>>>> Then in PCApply_YourShell(PC pc,Vec X,Vec Y) {
>>>>>>>> ...
>>>>>>>> PCGetOperators(pc,&A,PETSC_NULL,PETSC_NULL);
>>>>>>>> PCApply(ctx->fenpc,X,Y); /* apply first ordering */
>>>>>>>> VecScale(Y,-1.0);
>>>>>>>> MatMultAdd(A,Y,X,ctx->work); /* Compute fresh residual */
>>>>>>>> PCApply(ctx->fwnpc,ctx->work,Y);
>>>>>>>> }
>>>>>>>>> 
>>>>>>>>> I don't know if this is what you intended... is there a way to make
>>>>>>>>> KSP recompute the factorization with the different ordering without
>>>>>>>>> destroying and recreating the KSP?
>>>>>>>> 
>>>>>>>> Store both as in the code above.
>>>>>>>> 
>>>>>>>>> 
>>>>>>>>> Concerning Barry's option, js it possible to use the PETSc ilu
>>>>>>>>> factorization function (and eventually the ones provided by external
>>>>>>>>> libraries) inside my PCSHELL, passing it the various ordering I will
>>>>>>>>> need? If so, what's the best way to access it?
>>>>>>>> 
>>>>>>>> Does the code above answer this question?
>>>>>>> 
>>>>>>> 
>>>>>>> 
>>>>>>> --
>>>>>>> "[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
>>>> 
>>>> 
>>> 
>>> 
>>> 
>>> -- 
>>> "[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
>>