<div dir="ltr"><div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Mon, Oct 8, 2018 at 8:41 PM Ellen M. Price <<a href="mailto:ellen.price@cfa.harvard.edu">ellen.price@cfa.harvard.edu</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hi Matt,<br>
<br>
Thanks for the speedy reply. I do have a few more questions/elaborations<br>
on my previous correspondence.<br>
<br>
As for the source term (and the boundary condition, actually), both are<br>
only *known* at mesh points, but are assumed continuous everywhere. How<br>
does this change things? Basically, I want to assume a linear<br>
interpolation (I'm using linear elements) of the source term over each<br>
tetrahedral cell and then the integral should be some linear combination<br>
of mesh point values, with coefficients given by the inner product of<br>
basis functions? </blockquote><div><br></div><div>Its a little bit strange to give a forcing term this way. It means that you do not know</div><div>your problem independent of your discretization. For example, ex12 can use quadratic</div><div>elements, or cubic, or DG, so it does not assume any discretization when specifying the</div><div>problem.</div><div><br></div><div>You can do what you want, namely just pass in a rhs vector. However, that also means you</div><div>need to understand the mesh ordering, dof ordering, etc. Its not a great way to proceed.</div><div><br></div><div>How do you create these "mesh point" values in the first place?</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">The boundary term is actually given by a very nasty<br>
integral equation, which is probably the most expensive part of the<br>
simulation. So I only want to evaluate it at mesh points if possible.<br></blockquote><div><br></div><div>In order to maintain accuracy, you should really evaluate it at quadrature points (you can always</div><div>choose a low-order quadrature if you don't care). You can look at the Neumann term I use</div><div><br></div><div> <a href="https://bitbucket.org/petsc/petsc/src/760cff7d6b13b58cab0db4b57ca0fc45c2034c72/src/snes/examples/tutorials/ex12.c#lines-159">https://bitbucket.org/petsc/petsc/src/760cff7d6b13b58cab0db4b57ca0fc45c2034c72/src/snes/examples/tutorials/ex12.c#lines-159</a></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
You are correct, I forgot the \cdot n in that equation. Oops!<br>
<br>
I need to think more about the Jacobian and residual before I can really<br>
ask an intelligent question about them. I suppose that in practice, I<br>
need to figure out exactly which Jacobian terms I need to include, and<br>
which can be NULL. Some guidance on that would be immensely helpful.<br></blockquote><div><br></div><div>Put them all in. Always check the full physics and then use approximations as preconditioners.</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
If it helps, the equation I'm solving is just the Newtonian<br>
gravitational potential equation (on a nontrivial mesh and with a<br>
nontrivial, non-analytic mass density function).<br></blockquote><div><br></div><div>Hmm, FEM is not great for such rough source terms. At least not without some kind of nice adaptivity.</div><div><br></div><div>Also, I assume you are doing a boundary integral for boundary conditions. This is quickly going to</div><div>dominate the computational cost. I would figure out my strategy for that pretty early.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Ellen<br>
<br>
<br>
On 10/08/2018 06:28 PM, Matthew Knepley wrote:<br>
> On Mon, Oct 8, 2018 at 5:50 PM Ellen M. Price<br>
> <<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a> <mailto:<a href="mailto:ellen.price@cfa.harvard.edu" target="_blank">ellen.price@cfa.harvard.edu</a>>> wrote:<br>
> <br>
> Can anyone shed some light on SNES example 12?<br>
> <br>
> <a href="https://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex12.c.html" rel="noreferrer" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex12.c.html</a><br>
> <br>
> I am solving a Poisson problem, but with a source term defined on mesh<br>
> points, so I don't have an analytic form.<br>
> <br>
> <br>
> Can you explain more about this source term? It sounds like a bunch of<br>
> delta functions. That<br>
> would still work in this framework, but the convergence rate for a rhs<br>
> with these singularities<br>
> is reduced (this is a generic feature of FEM).<br>
> <br>
> <br>
> The part I am confused about<br>
> is the role of the functions f0, f1, and g3 that go into the DS.<br>
> <br>
> Suppose I have the Poisson problem nabla^2 u = f. Then the weak form is,<br>
> for test function phi,<br>
> <br>
> \int_{\partial \Omega} \phi \nabla u - \int_\Omega \nabla u \cdot \nabla<br>
> \phi = \int_\Omega \phi f<br>
> <br>
> <br>
> I think you need \int_{\partial \Omega} \phi \nabla u \cdot n instead.<br>
> <br>
> <br>
> Based on the documentation of PetscDSSetResidual(), I'm a little<br>
> confused about f0 and f1:<br>
> <br>
> \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec<br>
> f}_1(u, u_t, \nabla u, x, t)<br>
> <br>
> So f0 is, for the Poisson problem, zero, and f1 is just -grad u?<br>
> <br>
> <br>
> No. f1 is indeed -grad u, but f0 is -f.<br>
> <br>
> Note that ex12 solves -nabla^2 u = f since that operator is positive<br>
> definite.<br>
> <br>
> <br>
> Do we<br>
> need to do anything with the surface term, or is that handled<br>
> internally?<br>
> <br>
> <br>
> The surface term can be included using<br>
> <br>
> <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DT/PetscDSSetBdResidual.html" rel="noreferrer" target="_blank">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DT/PetscDSSetBdResidual.html</a><br>
> <br>
> and we use that in ex12 for inhomogeneous Neumann conditions.<br>
> <br>
> <br>
> Because neither of these terms are a surface integral term,<br>
> both are volume integrals.<br>
> <br>
> I can't figure out the purpose of g0, g1, g2, g3 in<br>
> PetscDSSetJacobian(), either. Is g supposed to be the source term?<br>
> <br>
> <br>
> No, it is the coefficient of the different combinations of basis<br>
> functions and their gradients.<br>
> For instance, the g3 term (grad phi and grad psi), g3 plays the same<br>
> role as the elasticity<br>
> tensor C in mechanics.<br>
> <br>
> <br>
> Is<br>
> psi a second basis function?<br>
> <br>
> <br>
> Yes. If the method is Galerkin, it comes from the same space as phi.<br>
> <br>
> <br>
> What even *is* the Jacobian in a FEM<br>
> problem?<br>
> <br>
> <br>
> The Frechet derivative of the residual.<br>
> <br>
> <br>
> I've only seen them formulated as linear problems in tutorials,<br>
> so I'm not even clear why we need a nonlinear solver or Jacobian...<br>
> <br>
> <br>
> The Frechet derivative of a linear operator is itself. ex12 also<br>
> incorporates<br>
> nonlinear variants of Poisson.<br>
> <br>
> <br>
> Any help is appreciated. I like the idea of a general framework for FEM,<br>
> as I think it will reduce human error on my part, but I also want to<br>
> understand what I'm doing.<br>
> <br>
> <br>
> If you have any more questions, feel free to keep mailing.<br>
> <br>
> Thanks,<br>
> <br>
> Matt<br>
> <br>
> <br>
> Ellen Price<br>
> <br>
> <br>
> <br>
> -- <br>
> What most experimenters take for granted before they begin their<br>
> experiments is infinitely more interesting than any results to which<br>
> their experiments lead.<br>
> -- Norbert Wiener<br>
> <br>
> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> <<a href="http://www.cse.buffalo.edu/%7Eknepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/%7Eknepley/</a>><br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div></div>