Yingjie Wu via petsc-users <[email protected]> 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.
