On Mon, Oct 18, 2021 at 3:27 AM Pierre Seize <[email protected]> wrote:
> On 15/10/21 22:30, Matthew Knepley wrote: > > On Fri, Oct 15, 2021 at 10:16 AM Pierre Seize <[email protected]> > wrote: > >> I read everything again, I think I did not understand you at first. The >> first solution is to modify the DAG, so that the rightmost cell is linked >> to the leftmost face, right ? To do that, do I have to manually edit the >> DAG (the mesh is read from a file) ? >> > Yes, the DAG would be modified if you want it for some particular mesh > that we cannot read automatically. For example, we can read periodic GMsh > meshes. > >> If so, the mesh connectivity is like the one of a torus, then how does it >> work with the cells/faces coordinates ? >> > You let the coordinate field be in a DG space, so that it can have jumps. > > I'm not sure I fully understand this, what I'll do is that I will > experiment with a periodic GMSH mesh. > We normally assume that the coordinate field on a mesh is continuous, which is why we associate it with vertices. However you could define the field on cells, representing it by discontinuous polynomials of degree 1. > Now I think the second method may be more straightforward. What's the idea >> ? Get the mapping with DMGetLocalToGlobalMapping, then create the mapping >> corresponding to the periodicity with ISLocalToGlobalMappingCreate, and >> finally ISLocalToGlobalMappingConcatenate ? I'm not sure this is the way, >> and I did not find something like DMSetLocalToGlobalMapping to restore the >> modified mapping. >> > It is more complicated. We make the LocalToGlobalMap by looking at the > PetscSection (essentially if gives function space information) and > deciding which unknowns are removed from the global space. > You would need to decide that unknowns constrained by periodicity are not > present in the global space. Actually, this is not hard. You just mark them > as constrained in the PetscSection, and all the layout > functions will function correctly. However, then the LocalToGlobalMap will > not be exactly right because the constrained unknowns will not be filled in > (just like Dirichlet conditions). You would augment the > map so that it fills those in by looking up their periodic counterparts. > Jed has argued for this type of periodicity. > > I also don't understand. One one side of my mesh I have : | ghost1 | > cell1 | cell2 | ... and on the other | cell_n-1 | cell_n | ghost_n |. Are > not the ghosts (from DMPlexConstructGhostCells) already constrained ? > I experimented on that too, I did: > > DMGetSectionSF > PetscSFGetGraph > augment the graph by adding the local ghost cells to the leaves, and the > correct remote "true" cells to the roots > PetscSFSetGraph > > and it seems to do what I want. Is this what you meant ? Is this a correct > way to use the PETSc objects ? Or is this just hacky and I'm lucky it works > ? > Yes, this will do what you want for field values, but not for coordinates. This is exactly what you would get if you just connected the topology. Thanks, Matt > Pierre > > > To me, the first kind is much more straightforward, but maybe this is > because I find the topology code more clear. > > Thanks, > > Matt > >> Pierre >> >> On 15/10/21 15:33, Pierre Seize wrote: >> >> When I first tried to handle the periodicity, I found the >> DMPlexCreateBoxMesh function (I cannot find the cylinder one). >> >> From reading the sources, I understand that we do some work either in >> DMPlexCreateCubeMesh_Internal or with DMSetPeriodicity. >> >> I tried to use DMSetPeriodicity before, for example with a 2x2 box on >> length 10. I did something like: >> const PetscReal maxCell[] = {2, 2}; >> const PetscReal L[] = {10, 10}; >> const DMBoundaryType bd[] = {DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC}; >> DMSetPeriodicity(dm, PETSC_TRUE, maxCell, L, bd); >> // or: >> DMSetPeriodicity(dm, PETSC_TRUE, NULL, L, bd); >> >> but it did not work: >> VecSet(X, 1); >> DMGetLocalVector(dm, &locX); >> VecZeroEntries(locX); >> DMGlobalToLocalBegin(dm, X, INSERT_VALUES, locX); >> DMGlobalToLocalEnd(dm, X, INSERT_VALUES, locX); >> VecView(locX, PETSC_VIEWER_STDOUT_WORLD); >> >> but the ghost cells values are all 0, only the real cells are 1. So I >> guess DMSetPeriodicity alone is not sufficient to handle the >> periodicity. Is there a way to do what I want ? That is set up my DMPlex in >> a way that DMGlobalToLocalBegin/DMGlobalToLocalEnd do exchange values >> between procs AND exchange the periodic values? >> >> >> Thanks for the help >> >> >> Pierre >> >> On 15/10/21 14:03, Matthew Knepley wrote: >> >> On Fri, Oct 15, 2021 at 7:31 AM Pierre Seize <[email protected]> >> wrote: >> >>> It makes sense, thank you. In fact, both ways seems better than my way. >>> The first one looks the most straightforward. Unfortunately I do not know >>> how to implement either of them. Could you please direct me to the >>> corresponding PETSc functions ? >>> >> The first way is implemented for example in DMPlexCreateBoxMesh() and >> DMPlexCreateCylinderMesh(). The second is not implemented since >> there did not seem to be a general way to do it. I would help if you >> wanted to try coding it up. >> >> Thanks, >> >> Matt >> >>> Pierre >>> >>> On 15/10/21 13:25, Matthew Knepley wrote: >>> >>> On Fri, Oct 15, 2021 at 7:08 AM Pierre Seize <[email protected]> >>> wrote: >>> >>>> Hi, >>>> >>>> I'm writing a code using PETSc to solve NS equations with FV on an >>>> unstructured mesh. Therefore I use DMPlex. >>>> >>>> Regarding periodicity, I manage to implement it this way: >>>> >>>> - for each couple of boundaries that is linked with periodicity, I >>>> create a buffer vector with an ISLocalToGlobalMapping >>>> >>>> - then, when I need to fill the ghost cells corresponding to the >>>> periodicity, the i "true" cell of the local vector fills the buffer >>>> vector on location i with VecSetValuesBlockedLocal, then >>>> VecAssemblyBegin/VecAssemblyEnd ensure each value is send to the >>>> correct >>>> location thanks to the mapping, then the i "ghost" cell of the local >>>> vector reads the vector on location i to get it's value. >>>> >>>> >>>> It works, but it seems to me there is a better way, with maybe PetscSF, >>>> VecScatter, or something I don't know yet. Does anyone have any advice ? >>>> >>> >>> There are at least two other ways to handle this. First, the method that >>> is advocated in >>> Plex is to actually make a periodic geometry, meaning connect the cells >>> that are meant >>> to be connected. Then, if you partition with overlap = 1, >>> PetscGlobalToLocal() will fill in >>> these cell values automatically. >>> >>> Second, you could use a non-periodic geometry, but alter the >>> LocalToGlobal map such >>> that the cells gets filled in anyway. Many codes use this scheme and it >>> is straightforward >>> with Plex just by augmenting the map it makes automatically. >>> >>> Does this make sense? >>> >>> Thanks, >>> >>> Matt >>> >>> >>>> Pierre Seize >>>> >>> -- >>> 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/%7Eknepley/> >>> >>> >>> >> >> -- >> 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/%7Eknepley/> >> >> >> >> > > -- > 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/%7Eknepley/> > > > -- 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/>
