<div dir="ltr"><div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Mon, Oct 8, 2018 at 5:50 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">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.</blockquote><div><br></div><div>Can you explain more about this source term? It sounds like a bunch of delta functions. That</div><div>would still work in this framework, but the convergence rate for a rhs with these singularities</div><div>is reduced (this is a generic feature of FEM).</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 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></blockquote><div><br></div><div>I think you need \int_{\partial \Omega} \phi \nabla u \cdot n instead.</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">
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?</blockquote><div><br></div><div>No. f1 is indeed -grad u, but f0 is -f.</div><div><br></div><div>Note that ex12 solves -nabla^2 u = f since that operator is positive definite.</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"> Do we<br>
need to do anything with the surface term, or is that handled<br>
internally?</blockquote><div><br></div><div>The surface term can be included using</div><div><br></div><div>  <a href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DT/PetscDSSetBdResidual.html">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DT/PetscDSSetBdResidual.html</a></div><div><br></div><div>and we use that in ex12 for inhomogeneous Neumann conditions.</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"> 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?</blockquote><div><br></div><div>No, it is the coefficient of the different combinations of basis functions and their gradients.</div><div>For instance, the g3 term (grad phi and grad psi), g3 plays the same role as the elasticity</div><div>tensor C in mechanics.</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"> Is<br>
psi a second basis function?</blockquote><div><br></div><div>Yes. If the method is Galerkin, it comes from the same space as phi.</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"> What even *is* the Jacobian in a FEM<br>
problem?</blockquote><div><br></div><div>The Frechet derivative of the residual.</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"> 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></blockquote><div><br></div><div>The Frechet derivative of a linear operator is itself. ex12 also incorporates</div><div>nonlinear variants of Poisson.</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">
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></blockquote><div><br></div><div>If you have any more questions, feel free to keep mailing.</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 Price<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>