[petsc-dev] integer divide by zero when using a two triangle DMPlex

[Matthew Knepley]

On Thu, Nov 21, 2013 at 7:58 PM, Geoffrey Irving <irving at naml.us> wrote:

> On Thu, Nov 21, 2013 at 4:51 PM, Matthew Knepley <knepley at gmail.com>
> wrote:
> > On Thu, Nov 21, 2013 at 6:30 PM, Geoffrey Irving <irving at naml.us> wrote:
> >>
> >> On Thu, Nov 21, 2013 at 4:11 PM, Matthew Knepley <knepley at gmail.com>
> >> wrote:
> >> > On Thu, Nov 21, 2013 at 1:48 PM, Geoffrey Irving <irving at naml.us>
> wrote:
> >> >>
> >> >> I'm getting the following crash trying to use a two triangle mesh
> with
> >> >> DMPlex FEM.  The problem is that the PetscFEGetQuadrature call on
> line
> >> >> 546 of DMPlexComputeResidualFEM produces a quadrature object with
> >> >> q.numPoints = 0, which we then use for division.
> >> >
> >> > Yes, it causes batchSize = 0. I think we could error is q.numPoints ==
> >> > 0.
> >> >
> >> >> I'm using the default elements, which based on the "fe dofs" printout
> >> >> are piecewise constant.
> >> >
> >> > You are supposed to create the quadrature (like ex12.c:382). The
> default
> >> > thing is completely empty.
> >>
> >> The code initializes quadrature just like ex12:
> >>
> >>   /* Create quadrature (with specified order if given) */
> >>   PetscQuadrature q;
> >>
> CHECK(PetscDTGaussJacobiQuadrature(dim,qorder>0?qorder:order,-1,1,&q));
> >>   CHECK(PetscFESetQuadrature(fe,q));
> >>
> >> It works fine with "-petscspace_order 1" or "2", and computes the
> >> right answer to machine precision for 2 since the exact solution is
> >> quadratic.  The problem is that if you ask for zeroth order
> >> quadrature, you get no quadrature points.  Note that I'm using the
> >> same order quadrature as the polynomial space (I copied the code from
> >> ex12).
> >
> > It appears that my use of the term 'order' is the problem. Jed does not
> like it either.
> > This is really the number of nodes.
>
> If you're referring to the "order" parameter to
> PetscDTGaussJacobiQuadrature, there's no sense in which it could have
> been meant to mean number of nodes.  The first line is
>

Yes, total number of nodes, accounting for dimension.


>     PetscInt       npoints = dim > 1 ? dim > 2 ? order*PetscSqr(order)
> : PetscSqr(order) : order;
>
> >> So arguably user error, but it would be nice if the defaults didn't
> >> cause integer division exceptions.
> >
> > Yes, we should have a check for numPoints > 0 in the integration
> function.
> >
> >> Is silently producing no quadrature points the right thing for
> >> PetscFESetQuadrature to do?  If so, can I check somewhere else to
> >> avoid the zero division?
> >
> > I do not know. I tried to allow this corner case, but I am not sure its
> > useful.
>
> I can't think of a reason someone would want a quadrature rule with no
> points, other than to indicate that the quadrature rule simply isn't
> defined.  I suppose that's a valid reason, though.
>
> So, I started adding checks to the integration routines, and it seemed
> pretty ridiculous (PetscFEGetQuadrature is called a lot).  It seems
> like the check should go either in PetscDTGaussJacobiQuadrature and
> friends or PetscFESetQuadrature.  Or maybe PetscFEGetQuadrature?  I'm
> not sure what the standards are here, or whether we want to support
> setting and detecting explicitly uninitialized quadrature rules.


A check would definitely go in PetscDTGaussJacobiQuadrature() if we
want to disallow empty quadrature, and I just not sure we want to do that,
since it makes things hard in parallel when you have empty processes. I
think I would rather put the high level check in IntegrateResidual().

   Matt


> Geoffrey
>



-- 
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20131122/2ccf2195/attachment.html>