The libMesh/MOOSE specific code that identifies dof indices for ISCreateGeneral is in DMooseGetEmbedding_Private. I can share that function (it's quite long) or more details if that could be helpful.
On Mon, Nov 7, 2022 at 10:55 AM Alexander Lindsay <[email protected]> wrote: > 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 <[email protected]> wrote: > >> On Thu, Nov 3, 2022 at 7:52 PM Alexander Lindsay < >> [email protected]> 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/> >> >
