[petsc-users] ILU and Block Gauss Seidel smoothing

[Gianluca Meneghello]

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).

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?

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