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

[Smith, Barry F.]

> On Nov 4, 2018, at 8:45 AM, Yingjie Wu via petsc-users <petsc-users at mcs.anl.gov> 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