[petsc-users] -snes_type test shows differences
(Rebecca) Xuefei YUAN
xy2102 at columbia.edu
Wed Feb 10 18:04:05 CST 2010
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
More information about the petsc-users
mailing list