<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Mon, Jun 5, 2017 at 12:18 PM, Lawrence Mitchell <span dir="ltr"><<a href="mailto:lawrence.mitchell@imperial.ac.uk" target="_blank">lawrence.mitchell@imperial.ac.uk</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
> On 2 Jun 2017, at 05:09, Matthew Knepley <<a href="mailto:knepley@gmail.com">knepley@gmail.com</a>> wrote:<br>
><br>
><br>
> Coming back to this, I think I understand the problem a little better.<br>
><br>
> Consider this mesh:<br>
><br>
> +----+<br>
> |\ 3 |<br>
> | \ |<br>
> |2 \ |<br>
> | \|<br>
> +----+<br>
> |\ 1 |<br>
> | \ |<br>
> |0 \ |<br>
> | \|<br>
> +----+<br>
><br>
> Let's say I run on 3 processes and the initial (non-overlapped) cell<br>
> partition is:<br>
><br>
> rank 0: cell 0<br>
> rank 1: cell 1 & 2<br>
> rank 2: cell 3<br>
><br>
> Now I'd like to grow the overlap such that any cell I can see through<br>
> a facet (and its closure) lives in the overlap.<br>
><br>
> So great, I just need a new adjacency relation that gathers<br>
> closure(support(point))<br>
><br>
> But, that isn't right, because now on rank 0, I will get a mesh that<br>
> looks like:<br>
><br>
> I do not understand why you think its not right. Toby and I are trying to push a formalism for<br>
> this understanding, in <a href="https://arxiv.org/abs/1508.02470" rel="noreferrer" target="_blank">https://arxiv.org/abs/1508.<wbr>02470</a>. So you say that if sigma is a dual<br>
> basis function associated with point p, then the support of its matching psi, sigma(psi) = 1<br>
> in the biorthogonal bases, is exactly star(p).<br>
><br>
> So, if you have no sigma sitting on your vertex, who cares if you put that extra edge and node<br>
> in. It will not affect the communication pattern for dofs. If you do, then shouldn't you be including<br>
> that edge?<br>
<br>
Hmm, I think we are talking at cross-purposes. Let me try and explain again where I am coming from:<br>
<br>
To do a FEM integral on some cell c I need:<br>
<br>
i) to evaluate coefficients at quadrature points (on the cell)<br>
ii) to evaluate basis functions at quadrature points (on the cell)<br>
<br>
for (i), I need all the dofs in closure(c).<br>
for (ii), I just need the definition of the finite element.<br>
<br>
To do a FEM integral on a facet f I need:<br>
<br>
i) to evaluate coefficients at quadrature points (on the facet)<br>
ii) to evaluate basis functions at quadrature points (on the facet)<br>
<br>
for (i), I need all the dofs in closure(support(f)).<br></blockquote><div><br></div><div>So this is a jump term, since I need both sides. Is it nonlinear? If its linear, things are easy</div><div>and just compute from both sides and add, but of course that does not work for nonlinear things.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
for (ii), I just need the definition of the finite element.<br>
<br>
So now, my model for how I want to global assembly of a facet integral is:<br>
<br>
loop over all facets:<br>
gather from global coefficient to local data<br>
evaluate coefficient at quad points<br>
perform integral<br>
local to global<br>
<br>
In parallel, I just make a partition of the facets (so that each facet is integrated exactly once).<br>
<br>
OK, so what data do I need in parallel?<br>
<br>
Exactly the dofs that correspond to closure(support(facet)) for all owned facets in my partition.<br>
<br>
So I was hoping to be able to grow a distributed cell partition by exactly that means: add in those remote cells which are in the support of owned facets (I'm happy if this is symmetrically grown, although I think I can do it with one-sided growth).<br>
<br>
So that's my rationale for wanting this "strange" adjacency. I can get enough adjacency by using the current FEM adjacency and filtering which entities I iterate over, but it seems a bit wasteful.<br></blockquote><div><br></div><div>I see that you have edges you do not necessarily want, but do they mess up your loop? It seems like you will not encounter them looping over facets.</div><div>This is exactly what happens to me in FV, where I just ignore them.</div><div><br></div><div>If you do really want to prune them, then I guess overriding the DMPlexGetAdjacency() as you propose is probably the best way. I would</div><div>be willing to put it in. Please send me a reminder email since this week is pretty heinous for me.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
For the currently implemented adjacencies, however, the code definitely assumes in various places that the closure of the cells on a partition covers all the points. And doing, say, DMPlexSetAdjacencyUseCone/<wbr>Closure(PETSC_FALSE) results in meshes where that is not true.<br>
<br>
See the below:<br>
<br>
Lawrence<br>
<br>
$ mpiexec -n 3 ./bork<br>
[0] face 7 has support: 0<br>
[0] face 8 has support: 0<br>
[0] face 9 has support: 0 1<br>
[0] face 10 has support: 1<br>
[0] face 11 has support: 1<br>
[0] face 12 has support:<br>
[1] face 10 has support: 0 2<br>
[1] face 11 has support: 0 1<br>
[1] face 12 has support: 0<br>
[1] face 13 has support: 1<br>
[1] face 14 has support: 2<br>
[1] face 15 has support: 2<br>
[1] face 16 has support: 1 3<br>
[1] face 17 has support: 3<br>
[1] face 18 has support: 3<br>
[2] face 7 has support: 0 1<br>
[2] face 8 has support: 0<br>
[2] face 9 has support: 0<br>
[2] face 10 has support: 1<br>
[2] face 11 has support:<br>
[2] face 12 has support: 1[2]PETSC ERROR: --------------------- Error Message ------------------------------<wbr>------------------------------<wbr>--<br>
[2]PETSC ERROR: Petsc has generated inconsistent data<br>
[2]PETSC ERROR: Number of depth 2 faces 6 does not match permuted nubmer 5<br>
[2]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>documentation/faq.html</a> for trouble shooting.<br>
[2]PETSC ERROR: Petsc Development GIT revision: v3.7.6-3857-gda66ab19e3 GIT Date: 2017-05-10 09:02:09 -0500<br>
[2]PETSC ERROR: ./bork on a arch-darwin-c-opt named yam-laptop.local by lmitche1 Mon Jun 5 18:15:31 2017<br>
[2]PETSC ERROR: Configure options --download-chaco=1 --download-ctetgen=1 --download-eigen --download-exodusii=1 --download-hdf5=1 --download-hypre=1 --download-metis=1 --download-ml=1 --download-mumps=1 --download-netcdf=1 --download-parmetis=1 --download-ptscotch=1 --download-scalapack=1 --download-triangle=1 --with-c2html=0 --with-debugging=0 --with-shared-libraries=1 PETSC_ARCH=arch-darwin-c-opt<br>
[2]PETSC ERROR: #1 DMPlexCreateOrderingClosure_<wbr>Static() line 41 in /Users/lmitche1/Documents/<wbr>work/src/deps/petsc/src/dm/<wbr>impls/plex/plexreorder.c<br>
<br>
[0]PETSC ERROR: --------------------- Error Message ------------------------------<wbr>------------------------------<wbr>--<br>
[0]PETSC ERROR: Petsc has generated inconsistent data<br>
[0]PETSC ERROR: Number of depth 2 faces 6 does not match permuted nubmer 5<br>
[0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>documentation/faq.html</a> for trouble shooting.<br>
[0]PETSC ERROR: Petsc Development GIT revision: v3.7.6-3857-gda66ab19e3 GIT Date: 2017-05-10 09:02:09 -0500<br>
[0]PETSC ERROR: ./bork on a arch-darwin-c-opt named yam-laptop.local by lmitche1 Mon Jun 5 18:15:31 2017<br>
[0]PETSC ERROR: Configure options --download-chaco=1 --download-ctetgen=1 --download-eigen --download-exodusii=1 --download-hdf5=1 --download-hypre=1 --download-metis=1 --download-ml=1 --download-mumps=1 --download-netcdf=1 --download-parmetis=1 --download-ptscotch=1 --download-scalapack=1 --download-triangle=1 --with-c2html=0 --with-debugging=0 --with-shared-libraries=1 PETSC_ARCH=arch-darwin-c-opt<br>
[0]PETSC ERROR: #1 DMPlexCreateOrderingClosure_<wbr>Static() line 41 in /Users/lmitche1/Documents/<wbr>work/src/deps/petsc/src/dm/<wbr>impls/plex/plexreorder.c<br>
[0]PETSC ERROR: #2 DMPlexGetOrdering() line 133 in /Users/lmitche1/Documents/<wbr>work/src/deps/petsc/src/dm/<wbr>impls/plex/plexreorder.c<br>
[0]PETSC ERROR: #3 main() line 87 in /Users/lmitche1/Documents/<wbr>work/src/petsc-doodles/bork.c<br>
[0]PETSC ERROR: No PETSc Option Table entries<br>
[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov-------<wbr>---<br>
[2]PETSC ERROR: #2 DMPlexGetOrdering() line 133 in /Users/lmitche1/Documents/<wbr>work/src/deps/petsc/src/dm/<wbr>impls/plex/plexreorder.c<br>
[2]PETSC ERROR: #3 main() line 87 in /Users/lmitche1/Documents/<wbr>work/src/petsc-doodles/bork.c<br>
[2]PETSC ERROR: No PETSc Option Table entries<br>
[2]PETSC ERROR: ----------------End of Error Message -------send entire error message to petsc-maint@mcs.anl.gov-------<wbr>---<br>
application called MPI_Abort(MPI_COMM_WORLD, 77) - process 0<br>
<br>
#include <petsc.h><br>
#include <petscsf.h><br>
<br>
int main(int argc, char **argv)<br>
{<br>
PetscErrorCode ierr;<br>
DM dm = NULL, dmParallel = NULL;<br>
PetscSF sf = NULL;<br>
PetscPartitioner partitioner = NULL;<br>
IS perm = NULL;<br>
const PetscReal coords[12] = {0, 0,<br>
0, 0.5,<br>
0, 1,<br>
1, 0,<br>
1, 0.5,<br>
1, 1};<br>
const PetscInt cells[12] = {0, 1, 3,<br>
1, 4, 3,<br>
1, 2, 4,<br>
2, 5, 4};<br>
<br>
const PetscInt sizes[3] = {1, 2, 1};<br>
const PetscInt points[4] = {0, 1, 2, 3};<br>
PetscInt fStart, fEnd;<br>
PetscMPIInt rank, size;<br>
MPI_Comm comm;<br>
<br>
PetscInitialize(&argc, &argv, NULL, NULL);<br>
<br>
comm = PETSC_COMM_WORLD;<br>
<br>
ierr = MPI_Comm_rank(comm, &rank); CHKERRQ(ierr);<br>
ierr = MPI_Comm_size(comm, &size); CHKERRQ(ierr);<br>
<br>
if (size != 3) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Requires 3 processes");<br>
<br>
if (!rank) {<br>
ierr = DMPlexCreateFromCellList(comm, 2, 4, 6, 3, PETSC_TRUE,<br>
cells, 2, coords, &dm); CHKERRQ(ierr);<br>
} else {<br>
ierr = DMPlexCreateFromCellList(comm, 2, 0, 0, 3, PETSC_TRUE,<br>
NULL, 2, NULL, &dm); CHKERRQ(ierr);<br>
}<br>
<br>
ierr = PetscPartitionerCreate(comm, &partitioner); CHKERRQ(ierr);<br>
ierr = PetscPartitionerSetType(<wbr>partitioner, PETSCPARTITIONERSHELL); CHKERRQ(ierr);<br>
if (!rank) {<br>
ierr = PetscPartitionerShellSetPartit<wbr>ion(partitioner, size, sizes, points); CHKERRQ(ierr);<br>
} else {<br>
ierr = PetscPartitionerShellSetPartit<wbr>ion(partitioner, 3, NULL, NULL); CHKERRQ(ierr);<br>
}<br>
ierr = DMPlexSetPartitioner(dm, partitioner); CHKERRQ(ierr);<br>
<br>
ierr = DMPlexDistribute(dm, 0, &sf, &dmParallel); CHKERRQ(ierr);<br>
ierr = DMDestroy(&dm); CHKERRQ(ierr);<br>
ierr = PetscSFDestroy(&sf); CHKERRQ(ierr);<br>
ierr = PetscPartitionerDestroy(&<wbr>partitioner); CHKERRQ(ierr);<br>
<br>
ierr = DMViewFromOptions(dmParallel, NULL, "-parallel_dm_view"); CHKERRQ(ierr);<br>
<br>
ierr = DMPlexSetAdjacencyUseCone(<wbr>dmParallel, PETSC_FALSE); CHKERRQ(ierr);<br>
ierr = DMPlexSetAdjacencyUseClosure(<wbr>dmParallel, PETSC_FALSE); CHKERRQ(ierr);<br>
<br>
ierr = DMPlexDistributeOverlap(<wbr>dmParallel, 1, &sf, &dm); CHKERRQ(ierr);<br>
<br>
ierr = DMDestroy(&dmParallel); CHKERRQ(ierr);<br>
ierr = PetscSFDestroy(&sf); CHKERRQ(ierr);<br>
<br>
ierr = DMViewFromOptions(dm, NULL, "-overlap_dm_view"); CHKERRQ(ierr);<br>
<br>
ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd); CHKERRQ(ierr);<br>
<br>
for (PetscInt f = fStart; f < fEnd; f++) {<br>
const PetscInt *support = NULL;<br>
PetscInt supportSize;<br>
ierr = DMPlexGetSupportSize(dm, f, &supportSize); CHKERRQ(ierr);<br>
ierr = DMPlexGetSupport(dm, f, &support); CHKERRQ(ierr);<br>
ierr = PetscSynchronizedPrintf(comm, "[%d] face %d has support:", rank, f); CHKERRQ(ierr);<br>
for (PetscInt p = 0; p < supportSize; p++) {<br>
ierr = PetscSynchronizedPrintf(comm, " %d", support[p]); CHKERRQ(ierr);<br>
}<br>
ierr = PetscSynchronizedPrintf(comm, "\n"); CHKERRQ(ierr);<br>
}<br>
ierr = PetscSynchronizedFlush(comm, PETSC_STDOUT); CHKERRQ(ierr);<br>
<br>
ierr = DMPlexGetOrdering(dm, MATORDERINGRCM, NULL, &perm); CHKERRQ(ierr);<br>
if (perm) {<br>
ierr = ISViewFromOptions(perm, NULL, "-ordering_is_view"); CHKERRQ(ierr);<br>
ierr = ISDestroy(&perm); CHKERRQ(ierr);<br>
}<br>
ierr = DMDestroy(&dm); CHKERRQ(ierr);<br>
ierr = PetscFinalize();<br>
<br>
return 0;<br>
}<br>
<br>
<br>
<br>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">http://www.caam.rice.edu/~mk51/</a><br></div></div></div>
</div></div>