[petsc-users] ILU and Block Gauss Seidel smoothing
Gianluca Meneghello
gianmail at gmail.com
Mon Jul 18 13:36:36 CDT 2011
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
>
More information about the petsc-users
mailing list