[petsc-users] Implementing FEM with PLEX/DS/SNES

[Matthew Knepley]

On Mon, Oct 8, 2018 at 8:41 PM Ellen M. Price <ellen.price at cfa.harvard.edu>
wrote:

> Hi Matt,
>
> Thanks for the speedy reply. I do have a few more questions/elaborations
> on my previous correspondence.
>
> As for the source term (and the boundary condition, actually), both are
> only *known* at mesh points, but are assumed continuous everywhere. How
> does this change things? Basically, I want to assume a linear
> interpolation (I'm using linear elements) of the source term over each
> tetrahedral cell and then the integral should be some linear combination
> of mesh point values, with coefficients given by the inner product of
> basis functions?


Its a little bit strange to give a forcing term this way. It means that you
do not know
your problem independent of your discretization. For example, ex12 can use
quadratic
elements, or cubic, or DG, so it does not assume any discretization when
specifying the
problem.

You can do what you want, namely just pass in a rhs vector. However, that
also means you
need to understand the mesh ordering, dof ordering, etc. Its not a great
way to proceed.

How do you create these "mesh point" values in the first place?


> The boundary term is actually given by a very nasty
> integral equation, which is probably the most expensive part of the
> simulation. So I only want to evaluate it at mesh points if possible.
>

In order to maintain accuracy, you should really evaluate it at quadrature
points (you can always
choose a low-order quadrature if you don't care). You can look at the
Neumann term I use


https://bitbucket.org/petsc/petsc/src/760cff7d6b13b58cab0db4b57ca0fc45c2034c72/src/snes/examples/tutorials/ex12.c#lines-159


> You are correct, I forgot the \cdot n in that equation. Oops!
>
> I need to think more about the Jacobian and residual before I can really
> ask an intelligent question about them. I suppose that in practice, I
> need to figure out exactly which Jacobian terms I need to include, and
> which can be NULL. Some guidance on that would be immensely helpful.
>

Put them all in. Always check the full physics and then use approximations
as preconditioners.


> If it helps, the equation I'm solving is just the Newtonian
> gravitational potential equation (on a nontrivial mesh and with a
> nontrivial, non-analytic mass density function).
>

Hmm, FEM is not great for such rough source terms. At least not without
some kind of nice adaptivity.

Also, I assume you are doing a boundary integral for boundary conditions.
This is quickly going to
dominate the computational cost. I would figure out my strategy for that
pretty early.

  Thanks,

     Matt


> Ellen
>
>
> On 10/08/2018 06:28 PM, Matthew Knepley wrote:
> > On Mon, Oct 8, 2018 at 5:50 PM Ellen M. Price
> > <ellen.price at cfa.harvard.edu <mailto:ellen.price at cfa.harvard.edu>>
> wrote:
> >
> >     Can anyone shed some light on SNES example 12?
> >
> >
> https://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex12.c.html
> >
> >     I am solving a Poisson problem, but with a source term defined on
> mesh
> >     points, so I don't have an analytic form.
> >
> >
> > Can you explain more about this source term? It sounds like a bunch of
> > delta functions. That
> > would still work in this framework, but the convergence rate for a rhs
> > with these singularities
> > is reduced (this is a generic feature of FEM).
> >
> >
> >     The part I am confused about
> >     is the role of the functions f0, f1, and g3 that go into the DS.
> >
> >     Suppose I have the Poisson problem nabla^2 u = f. Then the weak form
> is,
> >     for test function phi,
> >
> >     \int_{\partial \Omega} \phi \nabla u - \int_\Omega \nabla u \cdot
> \nabla
> >     \phi = \int_\Omega \phi f
> >
> >
> > I think you need \int_{\partial \Omega} \phi \nabla u \cdot n instead.
> >
> >
> >     Based on the documentation of PetscDSSetResidual(), I'm a little
> >     confused about f0 and f1:
> >
> >     \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec
> >     f}_1(u, u_t, \nabla u, x, t)
> >
> >     So f0 is, for the Poisson problem, zero, and f1 is just -grad u?
> >
> >
> > No. f1 is indeed -grad u, but f0 is -f.
> >
> > Note that ex12 solves -nabla^2 u = f since that operator is positive
> > definite.
> >
> >
> >     Do we
> >     need to do anything with the surface term, or is that handled
> >     internally?
> >
> >
> > The surface term can be included using
> >
> >
> https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DT/PetscDSSetBdResidual.html
> >
> > and we use that in ex12 for inhomogeneous Neumann conditions.
> >
> >
> >     Because neither of these terms are a surface integral term,
> >     both are volume integrals.
> >
> >     I can't figure out the purpose of g0, g1, g2, g3 in
> >     PetscDSSetJacobian(), either. Is g supposed to be the source term?
> >
> >
> > No, it is the coefficient of the different combinations of basis
> > functions and their gradients.
> > For instance, the g3 term (grad phi and grad psi), g3 plays the same
> > role as the elasticity
> > tensor C in mechanics.
> >
> >
> >     Is
> >     psi a second basis function?
> >
> >
> > Yes. If the method is Galerkin, it comes from the same space as phi.
> >
> >
> >     What even *is* the Jacobian in a FEM
> >     problem?
> >
> >
> > The Frechet derivative of the residual.
> >
> >
> >     I've only seen them formulated as linear problems in tutorials,
> >     so I'm not even clear why we need a nonlinear solver or Jacobian...
> >
> >
> > The Frechet derivative of a linear operator is itself. ex12 also
> > incorporates
> > nonlinear variants of Poisson.
> >
> >
> >     Any help is appreciated. I like the idea of a general framework for
> FEM,
> >     as I think it will reduce human error on my part, but I also want to
> >     understand what I'm doing.
> >
> >
> > If you have any more questions, feel free to keep mailing.
> >
> >   Thanks,
> >
> >      Matt
> >
> >
> >     Ellen Price
> >
> >
> >
> > --
> > 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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181008/5f45fd80/attachment.html>