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

[Ellen M. Price]

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? 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.

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.

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).

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/>