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
