[petsc-users] Problem with TS and PCMG

[Jed Brown]

Torquil Macdonald Sørensen <torquil at gmail.com> writes:

> If you mean the entire backtrace within GDB, I'm including it in this
> email. 

Perfect, thank you.  Cc'ing Barry and Peter who might also like to comment.

> I'm using PCMG. My minimal testcase program has everything
> hardcoded. I'm not providing any runtime options other than
> -start_in_debugger, so everything can be seen from the code. I've
> since made a couple of small changes that might influence line number
> statements from GDB, so I'm including the newest version in this
> message.
>
> Here is the GDB output of "bt full" and the corresponding program code:
>
> (gdb) bt full
> #1  0x00007ffdc080d172 in DMCoarsen (dm=0x1fa7df0, comm=0x7ffdbffdd9c0
> <ompi_mpi_comm_null>, dmc=0x2048f00) at dm.c:1811
>         ierr = 0
>         link = 0x7ffdc1af1640
>         __func__ = "DMCoarsen"
> #2  0x00007ffdc0c41c1e in PCSetUp_MG (pc=0x1ff7570) at mg.c:593
>         kdm = 0x1d5b370
>         p = 0xb453
>         rscale = 0x1ff7570
>         lu = PETSC_FALSE
>         svd = (PETSC_TRUE | unknown: 32764)
>         dB = 0x0
>         __func__ = "PCSetUp_MG"
>         tvec = 0x0
>         viewer = 0x0
>         mglevels = 0x2006980
>         ierr = 0
>         preonly = PETSC_FALSE
>         cholesky = (unknown: 3246112504)
>         mg = 0x20060c0
>         i = 0
>         n = 2
>         redundant = (unknown: 1211234)
>         use_amat = PETSC_TRUE
>         dms = 0x2048f00
>         cpc = 0x200f620
>         dump = PETSC_FALSE
>         opsset = PETSC_FALSE
>         dA = 0x7ffdc1af1640
>         uflag = DIFFERENT_NONZERO_PATTERN

Okay, this is the relevant bit of code.

  /* Skipping this for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs? */
  if (pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
    /* construct the interpolation from the DMs */
    Mat p;
    Vec rscale;
    ierr     = PetscMalloc1(n,&dms);CHKERRQ(ierr);
    dms[n-1] = pc->dm;
    for (i=n-2; i>-1; i--) {
      DMKSP kdm;
      ierr = DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);CHKERRQ(ierr);
      ierr = KSPSetDM(mglevels[i]->smoothd,dms[i]);CHKERRQ(ierr);
      if (mg->galerkin) {ierr = KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);CHKERRQ(ierr);}
      ierr = DMGetDMKSPWrite(dms[i],&kdm);CHKERRQ(ierr);
      /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
       * a bitwise OR of computing the matrix, RHS, and initial iterate. */
      kdm->ops->computerhs = NULL;
      kdm->rhsctx          = NULL;
      if (!mglevels[i+1]->interpolate) {
        ierr = DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);CHKERRQ(ierr);

I'm inclined to implement DMCoarsen_Shell() so that it carries through
to the levels.  It won't be used for interpolation because of the final
conditional (the user has already set mglevels[i+1]->interpolate).  This
is slightly annoying because it's one more function that users have to
implement.  I think the patch below will also work, though I'm nervous
about DMShellSetCreateMatrix().

diff --git i/src/dm/interface/dm.c w/src/dm/interface/dm.c
index ec60763..c987200 100644
--- i/src/dm/interface/dm.c
+++ w/src/dm/interface/dm.c
@@ -1943,7 +1943,11 @@ PetscErrorCode  DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
 
   PetscFunctionBegin;
   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
-  ierr                      = (*dm->ops->coarsen)(dm, comm, dmc);CHKERRQ(ierr);
+  if (dm->ops->coarsen) {
+    ierr = (*dm->ops->coarsen)(dm, comm, dmc);CHKERRQ(ierr);
+  } else {
+    ierr = DMShellCreate(comm,dmc);CHKERRQ(ierr);
+  }
   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
   ierr                      = PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);CHKERRQ(ierr);
   (*dmc)->ctx               = dm->ctx;

>     Mat restriction; // Multigrid restriction matrix, not the transpose
> of "interpolation"
>     ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nb_dof_coarse, nb_dof_fine,
> 1, 0, &restriction); CHKERRQ(ierr);
>     ierr = MatSetValue(restriction, 0, 0, 1.0, INSERT_VALUES);
> CHKERRQ(ierr);
>     ierr = MatSetValue(restriction, 1, 2, 1.0, INSERT_VALUES);
> CHKERRQ(ierr);
>     ierr = MatSetValue(restriction, 2, 4, 1.0, INSERT_VALUES);
> CHKERRQ(ierr);

This "restriction" has primary order 0 (it aliases high-frequency
harmonics without attenuation), and is thus not suitable for residual
restriction.  It is a great _state_ restriction operator, as are used by
FAS methods.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 835 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131218/eb8a1595/attachment.pgp>