[petsc-users] Local columns of A10 do not equal local rows of A00

[Alexander Lindsay]

I'm not sure exactly what you mean, but I'll try to give more details. We
have our own DM class (DM_Moose) and we set our own field and domain
decomposition routines:

  dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Moose;

  dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_Moose;


The field and domain decomposition routines are as follows (can see also at
https://github.com/idaholab/moose/blob/next/framework/src/utils/PetscDMMoose.C
):

static PetscErrorCode
DMCreateFieldDecomposition_Moose(
    DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
{
  PetscErrorCode ierr;
  DM_Moose * dmm = (DM_Moose *)(dm->data);

  PetscFunctionBegin;
  /* Only called after DMSetUp(). */
  if (!dmm->_splitlocs)
    PetscFunctionReturn(0);
  *len = dmm->_splitlocs->size();
  if (namelist)
  {
    ierr = PetscMalloc(*len * sizeof(char *), namelist);
    CHKERRQ(ierr);
  }
  if (islist)
  {
    ierr = PetscMalloc(*len * sizeof(IS), islist);
    CHKERRQ(ierr);
  }
  if (dmlist)
  {
    ierr = PetscMalloc(*len * sizeof(DM), dmlist);
    CHKERRQ(ierr);
  }
  for (const auto & dit : *(dmm->_splitlocs))
  {
    unsigned int d = dit.second;
    std::string dname = dit.first;
    DM_Moose::SplitInfo & dinfo = (*dmm->_splits)[dname];
    if (!dinfo._dm)
    {
      ierr = DMCreateMoose(((PetscObject)dm)->comm, *dmm->_nl, &dinfo._dm);
      CHKERRQ(ierr);
      ierr = PetscObjectSetOptionsPrefix((PetscObject)dinfo._dm,
((PetscObject)dm)->prefix);
      CHKERRQ(ierr);
      std::string suffix = std::string("fieldsplit_") + dname + "_";
      ierr = PetscObjectAppendOptionsPrefix((PetscObject)dinfo._dm,
suffix.c_str());
      CHKERRQ(ierr);
    }
    ierr = DMSetFromOptions(dinfo._dm);
    CHKERRQ(ierr);
    ierr = DMSetUp(dinfo._dm);
    CHKERRQ(ierr);
    if (namelist)
    {
      ierr = PetscStrallocpy(dname.c_str(), (*namelist) + d);
      CHKERRQ(ierr);
    }
    if (islist)
    {
      if (!dinfo._rembedding)
      {
        IS dembedding, lembedding;
        ierr = DMMooseGetEmbedding_Private(dinfo._dm, &dembedding);
        CHKERRQ(ierr);
        if (dmm->_embedding)
        {
          // Create a relative embedding into the parent's index space.
          ierr = ISEmbed(dembedding, dmm->_embedding, PETSC_TRUE,
&lembedding);
          CHKERRQ(ierr);
          const PetscInt * lindices;
          PetscInt len, dlen, llen, *rindices, off, i;
          ierr = ISGetLocalSize(dembedding, &dlen);
          CHKERRQ(ierr);
          ierr = ISGetLocalSize(lembedding, &llen);
          CHKERRQ(ierr);
          if (llen != dlen)
            SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to
embed split %D", d);
          ierr = ISDestroy(&dembedding);
          CHKERRQ(ierr);
          // Convert local embedding to global (but still relative)
embedding
          ierr = PetscMalloc(llen * sizeof(PetscInt), &rindices);
          CHKERRQ(ierr);
          ierr = ISGetIndices(lembedding, &lindices);
          CHKERRQ(ierr);
          ierr = PetscMemcpy(rindices, lindices, llen * sizeof(PetscInt));
          CHKERRQ(ierr);
          ierr = ISDestroy(&lembedding);
          CHKERRQ(ierr);
          // We could get the index offset from a corresponding global
vector, but subDMs don't yet
          // have global vectors
          ierr = ISGetLocalSize(dmm->_embedding, &len);
          CHKERRQ(ierr);

          ierr = MPI_Scan(&len,
                          &off,
                          1,
#ifdef PETSC_USE_64BIT_INDICES
                          MPI_LONG_LONG_INT,
#else
                          MPI_INT,
#endif
                          MPI_SUM,
                          ((PetscObject)dm)->comm);
          CHKERRQ(ierr);

          off -= len;
          for (i = 0; i < llen; ++i)
            rindices[i] += off;
          ierr = ISCreateGeneral(
              ((PetscObject)dm)->comm, llen, rindices, PETSC_OWN_POINTER,
&(dinfo._rembedding));
          CHKERRQ(ierr);
        }
        else
        {
          dinfo._rembedding = dembedding;
        }
      }
      ierr = PetscObjectReference((PetscObject)(dinfo._rembedding));
      CHKERRQ(ierr);
      (*islist)[d] = dinfo._rembedding;
    }
    if (dmlist)
    {
      ierr = PetscObjectReference((PetscObject)dinfo._dm);
      CHKERRQ(ierr);
      (*dmlist)[d] = dinfo._dm;
    }
  }
  PetscFunctionReturn(0);
}

static PetscErrorCode
DMCreateDomainDecomposition_Moose(
    DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS **
outerislist, DM ** dmlist)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  /* Use DMCreateFieldDecomposition_Moose() to obtain everything but
outerislist, which is currently
   * PETSC_NULL. */
  if (outerislist)
    *outerislist = PETSC_NULL; /* FIX: allow mesh-based overlap. */
  ierr = DMCreateFieldDecomposition_Moose(dm, len, namelist, innerislist,
dmlist);
  CHKERRQ(ierr);
  PetscFunctionReturn(0);
}



On Thu, Nov 3, 2022 at 5:19 PM Matthew Knepley <knepley at gmail.com> wrote:

> On Thu, Nov 3, 2022 at 7:52 PM Alexander Lindsay <alexlindsay239 at gmail.com>
> wrote:
>
>> I have errors on quite a few (but not all) processes of the like
>>
>> [1]PETSC ERROR: --------------------- Error Message
>> --------------------------------------------------------------
>> [1]PETSC ERROR: Nonconforming object sizes
>> [1]PETSC ERROR: Local columns of A10 4137 do not equal local rows of A00
>> 4129
>>
>> when performing field splits. We (MOOSE) have some code for identifying
>> the index sets for each split. However, the code was written by some
>> authors who are no longer with us. Normally I would chase this down in a
>> debugger, but this error only seems to crop up for pretty complex and large
>> meshes. If anyone has an idea for what we might be doing wrong, that might
>> help me chase this down faster. I guess intuitively I'm pretty perplexed
>> that we could get ourselves into this pickle as it almost appears that we
>> have two different local dof index counts for a given block (0 in this
>> case). More background, if helpful, can be found in
>> https://github.com/idaholab/moose/issues/22359 as well as
>> https://github.com/idaholab/moose/discussions/22468.
>>
>
> How are you specifying the blocks? I would not have thought this was
> possible.
>
>   Thanks,
>
>      Matt
>
>
>> I should note that we are currently running with 3.16.6 as our PETSc
>> submodule hash (we are talking about updating to 3.18 soon).
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> https://www.cse.buffalo.edu/~knepley/
> <http://www.cse.buffalo.edu/~knepley/>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20221107/fd3f842f/attachment.html>