[petsc-users] -snes_type test shows differences

[(Rebecca) Xuefei YUAN]

Dear Jed,

Yes, I paste everything from the code, and it is only the  
value[4]([j,i]) wrong in the Jacobian.

I tried several different initials for f, but I always had such  
differences at the same point.

I doubt my analytic term for value[4] first, but I double checked and  
still cannot figure out what is wrong.

Thanks very much for your kind reply!

Cheers,

Rebecca
Quoting Jed Brown <jed at 59A2.org>:

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



-- 
(Rebecca) Xuefei YUAN
Department of Applied Physics and Applied Mathematics
Columbia University
Tel:917-399-8032
www.columbia.edu/~xy2102