Actual source code: ex25.c
1: static const char help[] ="Minimum surface problem in 2D.\n\
2: Uses 2-dimensional distributed arrays.\n\
3: \n\
4: Solves the linear systems via multilevel methods \n\
5: \n\n";
7: /*T
8: Concepts: SNES^solving a system of nonlinear equations
9: Concepts: DMDA^using distributed arrays
10: Concepts: multigrid;
11: Processors: n
12: T*/
14: /*
16: This example models the partial differential equation
18: - Div((1 + ||GRAD T||^2)^(1/2) (GRAD T)) = 0.
20: in the unit square, which is uniformly discretized in each of x and
21: y in this simple encoding. The degrees of freedom are vertex centered
23: A finite difference approximation with the usual 5-point stencil
24: is used to discretize the boundary value problem to obtain a
25: nonlinear system of equations.
27: */
29: #include <petscsnes.h>
30: #include <petscdm.h>
31: #include <petscdmda.h>
33: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,void*);
35: int main(int argc,char **argv)
36: {
37: SNES snes;
38: PetscInt its,lits;
39: PetscReal litspit;
40: DM da;
42: PetscInitialize(&argc,&argv,NULL,help);
43: /*
44: Set the DMDA (grid structure) for the grids.
45: */
46: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
47: DMSetFromOptions(da);
48: DMSetUp(da);
49: DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,NULL);
50: SNESCreate(PETSC_COMM_WORLD,&snes);
51: SNESSetDM(snes,da);
52: DMDestroy(&da);
54: SNESSetFromOptions(snes);
56: SNESSolve(snes,0,0);
57: SNESGetIterationNumber(snes,&its);
58: SNESGetLinearSolveIterations(snes,&lits);
59: litspit = ((PetscReal)lits)/((PetscReal)its);
60: PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);
61: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
62: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit);
64: SNESDestroy(&snes);
65: PetscFinalize();
66: return 0;
67: }
69: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr)
70: {
71: PetscInt i,j;
72: PetscScalar hx,hy;
73: PetscScalar gradup,graddown,gradleft,gradright,gradx,grady;
74: PetscScalar coeffup,coeffdown,coeffleft,coeffright;
77: hx = 1.0/(PetscReal)(info->mx-1); hy = 1.0/(PetscReal)(info->my-1);
79: /* Evaluate function */
80: for (j=info->ys; j<info->ys+info->ym; j++) {
81: for (i=info->xs; i<info->xs+info->xm; i++) {
83: if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {
84: f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
85: } else {
87: gradup = (t[j+1][i] - t[j][i])/hy;
88: graddown = (t[j][i] - t[j-1][i])/hy;
89: gradright = (t[j][i+1] - t[j][i])/hx;
90: gradleft = (t[j][i] - t[j][i-1])/hx;
92: gradx = .5*(t[j][i+1] - t[j][i-1])/hx;
93: grady = .5*(t[j+1][i] - t[j-1][i])/hy;
95: coeffup = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
96: coeffdown = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);
98: coeffleft = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
99: coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);
101: f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
102: }
104: }
105: }
106: return 0;
107: }
109: /*TEST
111: test:
112: args: -pc_type mg -da_refine 1 -ksp_type fgmres
114: test:
115: suffix: 2
116: nsize: 2
117: args: -pc_type mg -da_refine 1 -ksp_type fgmres
119: test:
120: suffix: 3
121: nsize: 2
122: args: -pc_type mg -da_refine 1 -ksp_type fgmres -snes_type newtontrdc -snes_trdc_use_cauchy false
124: test:
125: suffix: 4
126: nsize: 2
127: args: -pc_type mg -da_refine 1 -ksp_type fgmres -snes_type newtontrdc
128: filter: sed -e "s/SNES iterations = 1[1-3]/SNES iterations = 13/g" |sed -e "s/Linear iterations = 2[8-9]/Linear iterations = 29/g" |sed -e "s/Linear iterations = 3[0-1]/Linear iterations = 29/g"
130: TEST*/