> On Nov 4, 2018, at 8:45 AM, Yingjie Wu via petsc-users
> <[email protected]> wrote:
>
> 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?
>
> 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?
NonlinearGS() does iterations of point-block nonlinear Gauss-Seidel (which
is used as the smoother for nonlinear multigrid). Thus at each grid point it
does Newton's method (this is the inner loop for (l = 0; l < max_its &&
!ptconverged; l++) {). In nonlinear Gauss-Seidel smoothing only the diagonal
block of the Jacobian is needed, the off diagonal block are never needed.
Barry
>
> Thanks for your continuous help,
> Yingjie
