[petsc-users] ILU and Block Gauss Seidel smoothing

[Jed Brown]

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?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20110701/67adb3be/attachment-0001.htm>