[petsc-users] Some Problems about SNES EX19.c

[Jed Brown]

Yingjie Wu via petsc-users <petsc-users at mcs.anl.gov> writes:

> Dear Petsc developer:
> Hi,
> Recently, I am very interested in the ex19 example in SNES, which uses NGS
> method to solve the non-linear equations, which may be the method I need to
> use in the future. I have some doubts about the program.
> 1.
>
> 82: typedef struct {
> 83:   PetscScalar u,v,omega,temp;
> 84: } Field
> 86: PetscErrorCode FormFunctionLocal(DMDALocalInfo*,Field**,Field**,void*);
> …
> 150:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode
> (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,&user);
>
>
> I looked at PETSc manualpage:
>
> For PetscErrorCode (*func)(DMDALocalInfo *info,void *x, void *f, void *ctx),
> info - DMDALocalInfo defining the subdomain to evaluate the residual on
> x - dimensional pointer to state at which to evaluate residual (e.g.
> PetscScalar *x or **x or ***x)
> f - dimensional pointer to residual, write the residual here (e.g.
> PetscScalar *f or **f or ***f)
> ctx - optional context passed above
>
> In the function FormFunctionLocal, the second and third parameters should
> be pointers to PetscScalar, where pointers are directly used to point to
> Field. Why can we use them here?            Although I know that there are
> four degrees of freedom in DM objects, how can I ensure that the program
> correctly corresponds to variables in Field?

DMDA "interlaces" the fields, so this memory layout is guaranteed.

> 2.
> In the NGS subroutine,
>
> 530: dfudu = 2.0* (hydhx + hxdhy);
>
> But in the residual function:
>
> 526: u = x[j][i].u;
> 527: uxx = (2.0*u - x[j][i-1].u - x[j][i+1].u) *hydhx;
> 528: uyy = (2.0*u - x[j-1][i].u - x[j+1][i].u) *hxdhy;
> 529: Fu = uxx + uyy -.5* (x[j+1][i].omega-x[j-1][i].omega) *hx - bjiu;
>
> / * invert the system:
> 572:              [ dfu / du     0        0        0    ][yu] = [fu]
> 573:              [     0    dfv / dv     0        0    ][yv]   [fv]
> 574:              [ dfo / du dfo / dv dfo / do     0    ][yo]   [fo]
> 575:              [ dft / du dft / dv     0    dft / dt ][yt]   [ft]
> 576:              by simple back-substitution
> 577: * /
>
> It is known that the residual function fu [j] [i] is a function of five
> variables (x [j] [i].u, x [j] [i-1].u, x [j] [i+1].u, x [j-1] [i].u, x
> [j+1] [i] [i].u) (it is same for analytic Jacobian matrix). But in this
> program, only the central grid is used to solve the partial derivatives.
> Why do we choose to do so?            In my understanding, the sub-matrix
> 'dfudu' is a five-diagonal matrix, but it is processed into a diagonal
> matrix in the program. Will this affect the accuracy of the solution?

This may be a misunderstanding.  The residual calculates a vector as a
nonlinear function of the input.  The Jacobian matrix is not assembled
in this example, but can be computed efficiently by finite differencing
using coloring.  Each row has 5*4=20 formal nonzeros and computing the
matrix by coloring requires somewhere around that number of residual
evaluations.