[petsc-users] node DG with DMPlex
Yann Jobic
yann.jobic at univ-amu.fr
Fri Mar 27 13:27:50 CDT 2020
Hi matt,
Many thanks for the help !!
Le 3/26/2020 à 11:38 PM, Matthew Knepley a écrit :
>
> 17 18
> 7-----8-----9-----14----15
> | | |
> 16,3 4 5,6 11 12,13
> | | |
> 1-----2-----3-----9-----10
> 19 20
>
> so each face gets 2 dofs, one for each cell.When doing a cell integral,
> you only use the dof that is for that cell.
I think i understood. I put the beginning of the code for dof management
in attachment. i defined :
numDof[0*(user.dim+1)+user.dim] = internalDof; /* internalDof
defined on cells */
numDof[0*(user.dim+1)+user.dim-1] = nbDofFaceEle; /* nbDofFaceEle
defined on faces */
This way, the dof faces are not duplicated. They are defined on the
face. (this part is commented in the source file)
Am i right ?
> The local-to-global would update the face dofs, so you would get each
side.
>
This way, a local-to-global will update the face dofs ?
> There is a reordering when you extract the closure. I have written one
> for spectral elements. We would need
> another here that ordered all the "other" face dofs to the end.
You mean using PetscSectionSetPermutation ?
Or rewrite DMPlexSetClosurePermutationTensor in order to put all the
face dofs at the end ?
I don't understand this part. I'm looking at the documentation for
DMPlexSetClosurePermutationTensor.
https://www.mcs.anl.gov/petsc/petsc-current/src/dm/impls/plex/plex.c.html#DMPlexSetClosurePermutationTensor
It says that
----------------------
The closure in BFS ordering works through height strata (cells, edges,
vertices) to produce the ordering
.vb
0 1 2 3 8 9 14 15 11 10 13 12 4 5 7 6
.ve
--------------------------
which is what we want no ?
>
> This seems a little complicated to me. Do you know how Andreas Klockner
> does it in Hedge? Or Tim Warburton?
I'm reading his book (nodal discontinuous galerkin methods)
> I just want to make sure I am not missing an elegant way to handle this.
>
> Here field 1 is synchronised with field 0, locally.
> But the external Face dof, field 2, have to be synchronised with the
> values of the adjacent cell.
> Is it possible to use something like DMPlexGetFaceFields ?
> Is there an example of such use of PetscSection and synchronisation
> process ?
>
> For the parallel part, should i use PetscSF object ?
>
>
> In parallel, integrals would be summed into the global vector, so each
> side has a 0 for the other face dof and the right contribution
> for its face dof. Then both sides get both solution dofs. It seems to
> work in my head.
>
> I read your article "Mesh Algorithms for PDE with Sieve I: Mesh
> Distribution". But it's refereeing to Matthew G. Knepley and Dmitry A.
> Karpeev. Sieve implementation.
> Technical Report ANL/MCS to appear, Argonne National Laboratory,
> January 2008.
> I couldn't find it. It is freely available ?
>
>
> Don't bother reading that. There are later ones:
>
> There are two pretty good sources:
>
> https://arxiv.org/abs/1505.04633
> https://arxiv.org/abs/1506.06194
>
> The last one is a follow-on to this paper
>
> https://arxiv.org/abs/0908.4427
>
> Thanks,
>
> Matt
>
> >
> > I don't think you need extra vertices, > or coordinates, and for
> output I
> > recommend using DMPlexProject() to get
> > the solution in some space that can be plotted like P1, or
> anything else
> > supported by your visualization.
>
> I would like to use DMplex as much as i can, as i would in the future
> refine locally the mesh.
>
> I hope you're good in this difficult situation (covid19),
>
> Best regards,
>
> Yann
>
> >
> > Thanks,
> >
> > Matt
> >
> > >
> > > We have an implementation of spectral element ordering
> > >
> >
> (https://gitlab.com/petsc/petsc/-/blob/master/src/dm/impls/plex/examples/tutorials/ex6.c).
> >
> > > Those share
> > > the whole element boundary.
> > >
> > > 2) What ghosts do you need?
> > In order to compute the numerical fluxes of one element, i
> need the
> > values of the surrounding nodes connected to the adjacent
> elements.
> > >
> > > 3) You want to store real space coordinates for a
> quadrature?
> > It should be basically the same as PetscFE of higher order.
> > I add some vertex needed to compute a polynomal solution of
> the desired
> > order. That means that if i have a N, order of the local
> approximation,
> > i need 0.5*(N+1)*(N+2) vertex to store in the DMPlex (in 2D), in
> > order to :
> > 1) have the correct number of dof
> > 2) use ghost nodes to sync the values of the
> vertex/edge/facet for
> > 1D/2D/3D problem
> > 2) save correctly the solution
> >
> > Does it make sense to you ?
> >
> > Maybe like
> >
> https://www.mcs.anl.gov/petsc/petsc-current/src/ts/examples/tutorials/ex11.c.html
> > With the use of the function SplitFaces, which i didn't fully
> > understood
> > so far.
> >
> > Thanks,
> >
> > Yann
> >
> > >
> > > We usually define a quadrature on the reference
> element once.
> > >
> > > Thanks,
> > >
> > > Matt
> > >
> > > I found elements of answers in those threads :
> > >
> >
> https://lists.mcs.anl.gov/pipermail/petsc-users/2016-August/029985.html
> > >
> >
> https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2019-October/039581.html
> > >
> > > However, it's not clear for me where to begin.
> > >
> > > Quoting Matt, i should :
> > > " DMGetCoordinateDM(dm, &cdm);
> > > <Set field information into cdm instead of dm>
> > > DMCreateLocalVector(cdm, &coordinatesLocal);
> > > <Fill in higher order coordinate values>
> > > DMSetCoordinatesLocal(dm, coordinatesLocal);"
> > >
> > > However, i will not create ghost nodes this way. And
> i'm not
> > sure to
> > > keep the good ordering.
> > > This part should be implemented in the PetscFE
> interface, for
> > high
> > > order
> > > discrete solutions.
> > > I did not succeed in finding the correct part of the
> source
> > doing it.
> > >
> > > Could you please give me some hint to begin correctly
> thoses
> > tasks ?
> > >
> > > Thanks,
> > >
> > > Yann
> > >
> > >
> > >
> > > --
> > > 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/>
> >
> >
> >
> > --
> > 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/>
>
>
>
> --
> 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 --------------
static char help[] = "dof management with dmplex for DG problems, for simplices \n\n";
#include <petscdmplex.h>
typedef struct {
PetscInt dim; /* Topological problem dimension */
PetscInt Nf; /* Number of fields */
PetscInt Nc[1]; /* Number of components per field */
PetscInt N; /* Order of polymomials used for approximation */
} AppCtx;
static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
PetscBool flg;
PetscErrorCode ierr;
PetscFunctionBeginUser;
options->dim = 2;
options->Nf = 1;
options->Nc[0] = 1;
options->N = 4;
ierr = PetscOptionsBegin(comm, "", "dof management with dmplex for DG problems (simplices) Options", "DMPLEX");CHKERRQ(ierr);
ierr = PetscOptionsRangeInt("-dim", "Problem dimension", "dof_dg.c", options->dim, &options->dim, NULL,1,3);CHKERRQ(ierr);
ierr = PetscOptionsInt("-order", "Order of polymomials used for approximation", "dof_dg.c", options->N, &options->N, &flg);CHKERRQ(ierr);
ierr = PetscOptionsEnd();
PetscFunctionReturn(0);
}
int main(int argc, char **argv) {
DM dm;
PetscSection s;
Vec u;
PetscViewer viewer;
AppCtx user;
PetscInt cells[3] = {2, 2, 2};
PetscReal lower[3] = {-1,-1,-1};
PetscReal upper[3] = {1,1,1};
/* DG specification of the PetscSection */
/* order of poly approx is stored in AppCtx user : user.N */
PetscInt dofPerEle=0; /* Number of dof per cell */
PetscInt internalDof; /* dof defined on the internal cell */
/* so without the faces dof */
PetscInt nbFace=0; /* Number of faces per cell */
PetscInt nbDofFaceEle=0; /* Number of dof per face */
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr;
ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, user.dim, PETSC_FALSE, cells, lower, upper, NULL, PETSC_TRUE, &dm);CHKERRQ(ierr);
ierr = DMSetFromOptions(dm);CHKERRQ(ierr);
ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr);
/* Total number of dof per element, number of faces per element and number of dof per face */
if (user.dim == 2) {
dofPerEle = (user.N+1)*(user.N+2)/2;
nbFace = 3;
nbDofFaceEle = user.N + 1;
internalDof = dofPerEle - nbFace*nbDofFaceEle + 3;
} else if (user.dim == 3) {
dofPerEle = (user.N+1)*(user.N+2)*(user.N+3)/6;
nbFace = 4;
nbDofFaceEle = (user.N+1)*(user.N+2)/2;
internalDof = dofPerEle - nbFace*nbDofFaceEle + (user.N + 1)*6 - 4;
}
PetscPrintf(MPI_COMM_WORLD,"\ndof per cell : %d\n",dofPerEle);
PetscPrintf(MPI_COMM_WORLD,"dof per face : %d\n",nbDofFaceEle);
PetscPrintf(MPI_COMM_WORLD,"internal dof : %d\n",internalDof);
{
PetscInt *numDof, d;
/* numDof[f*(dim+1)+d] gives the number of dof for field f on */
/* points of dimension d. For instance, numDof[1] is the number */
/* of dof for field 0 on each edge. */
ierr = PetscMalloc1(user.Nf*(user.dim+1), &numDof);CHKERRQ(ierr);
/* fields 0, written generally to add easily more fields */
for (d = 0; d < user.Nf*(user.dim+1); ++d) numDof[d] = 0;
/* we put only the interior dof for not duplicating the face nodes */
numDof[0*(user.dim+1)+user.dim] = internalDof; /* internalDof defined on cells */
numDof[0*(user.dim+1)+user.dim-1] = nbDofFaceEle; /* nbDofFaceEle defined on faces */
ierr = DMSetNumFields(dm, user.Nf);CHKERRQ(ierr);
ierr = DMPlexCreateSection(dm, NULL, user.Nc, numDof, 0, NULL, NULL, NULL, NULL, &s);CHKERRQ(ierr);
ierr = PetscFree(numDof);CHKERRQ(ierr);
}
/* Name the Field variable */
ierr = PetscSectionSetFieldName(s, 0, "u");CHKERRQ(ierr);
/* Tell the DM to use this data layout */
ierr = DMSetLocalSection(dm, s);CHKERRQ(ierr);
/* Create a Vec with this layout and view it */
ierr = DMGetGlobalVector(dm, &u);CHKERRQ(ierr);
ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
ierr = PetscViewerFileSetName(viewer, "sol.vtk");CHKERRQ(ierr);
ierr = VecView(u, viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = DMRestoreGlobalVector(dm, &u);CHKERRQ(ierr);
/* Cleanup */
ierr = PetscSectionDestroy(&s);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
More information about the petsc-users
mailing list