[petsc-users] -snes_type test shows differences

[Jed Brown]

First, sorry for my off-hand comment.  Looking at your code more
carefully, it looks like you have the correct Newton linearization.  Did
you paste everything from the source?  I assume the indexing is correct
since you say it's only the diagonal value that's wrong.

On Wed, 10 Feb 2010 15:48:57 -0500, "(Rebecca) Xuefei YUAN" <xy2102 at columbia.edu> wrote:
> Dear Jed,
> 
> What if my F has the form(the full version):
> 
> F(f) = fxx+fyy+fxx*fyy-fxy*fxy-g = 0,
> 
> how should my Newton iteration look like?
> 
> J(f)h = -F(f) = -(fxx+fyy+fxx*hyy+fyy*hxx-fxy*hxy-fxy*hxy-g)?
> 
> I know h is the update, but what are hxx,hyy?

Those are derivatives of the update (which specify where things end up
in your stencil).  I find this quite error-prone to do by hand [1], but
symbolic algebra can do it very well.  As a simple example, name your
points North, South, East, West, Center (N,S,E,W,C).  Then the
discretization for the original function looks like

  from sympy import *
  from sympy.abc import *                     # defines the 1-letter symbols
  NW,NE,SW,SE = symbols('NW','NE','SW','SE')  # define symbols for the corner points
  fxx = (E+W-2*C)/h**2
  fyy = (N+S-2*C)/h**2
  fxy = ((NW-NE)/(2*h) - (SW-SE)/(2*h))/(2*h)
  F = (fxx + fyy + fxx*fyy - fxy*fxy) * h**2  # Final h**2 is just scaling

Now you can get your Jacobian with

  [simplify(diff(F,sym)) for sym in [NW,N,NE,W,C,E,SW,S,SE]]

which produces

  [(NE + SW - NW - SE)/(8*h**2),
   (E + W - 2*C + h**2)/h**2,
   (NW + SE - NE - SW)/(8*h**2),
   (N + S - 2*C + h**2)/h**2,
   (-2*E - 2*N - 2*S - 2*W + 8*C - 4*h**2)/h**2,
   (N + S - 2*C + h**2)/h**2,
   (NW + SE - NE - SW)/(8*h**2),
   (E + W - 2*C + h**2)/h**2,
   (NE + SW - NW - SE)/(8*h**2)]


You can get sympy to actually write the C code by defining symbols as

  NW,N = symbols(['x[j+1][i-1].field','x[j+1][i].field']) # ...

and using sympy.ccode(diff(F,sym)) to format the output.  Since the
Python code is so short, I often just paste it into the C code as a
comment to make it easier to change the physics later.  If you don't
trust the compiler's CSE, you can get efficient representations using
sympy.cse (in principle, more transformations are possible in the
symbolic world).

  http://sympy.org

Of course Maple/Mathematica/Maxima can also do this, but sympy is
lightweight and free.

I hope this helps.


Jed

[1] Error-prone for finite-difference methods.  I think it's usually
worth doing by hand for Galerkin methods because it can be done entirely
in the continuum setting and gives you insight about the problem.