Actual source code: ex25.c
2: static char help[] ="Minimum surface problem\n\
3: Uses 2-dimensional distributed arrays.\n\
4: \n\
5: Solves the linear systems via multilevel methods \n\
6: \n\n";
8: /*T
9: Concepts: SNES^solving a system of nonlinear equations
10: Concepts: DA^using distributed arrays
11: Concepts: multigrid;
12: Processors: n
13: T*/
15: /*
16:
17: This example models the partial differential equation
18:
19: - Div((1 + ||GRAD T||^2)^(1/2) (GRAD T)) = 0.
20:
21:
22: in the unit square, which is uniformly discretized in each of x and
23: y in this simple encoding. The degrees of freedom are vertex centered
24:
25: A finite difference approximation with the usual 5-point stencil
26: is used to discretize the boundary value problem to obtain a
27: nonlinear system of equations.
28:
29: */
31: #include petscsnes.h
32: #include petscda.h
33: #include petscmg.h
34: #include petscdmmg.h
41: int main(int argc,char **argv)
42: {
43: DMMG *dmmg;
44: SNES snes;
46: PetscInt its,lits;
47: PetscReal litspit;
48: DA da;
50: PetscInitialize(&argc,&argv,PETSC_NULL,help);
53: /*
54: Create the multilevel DA data structure
55: */
56: DMMGCreate(PETSC_COMM_WORLD,3,0,&dmmg);
58: /*
59: Set the DA (grid structure) for the grids.
60: */
61: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,-5,-5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da);
62: DMMGSetDM(dmmg,(DM)da);
63: DADestroy(da);
65: /*
66: Process adiC(36): FormFunctionLocal FormFunctionLocali
68: Create the nonlinear solver, and tell the DMMG structure to use it
69: */
70: /* DMMGSetSNES(dmmg,FormFunction,0); */
71: DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,ad_FormFunctionLocal,0);
72: DMMGSetFromOptions(dmmg);
74: /*
75: PreLoadBegin() means that the following section of code is run twice. The first time
76: through the flag PreLoading is on this the nonlinear solver is only run for a single step.
77: The second time through (the actually timed code) the maximum iterations is set to 10
78: Preload of the executable is done to eliminate from the timing the time spent bring the
79: executable into memory from disk (paging in).
80: */
81: PreLoadBegin(PETSC_TRUE,"Solve");
82: DMMGSolve(dmmg);
83: PreLoadEnd();
84: snes = DMMGGetSNES(dmmg);
85: SNESGetIterationNumber(snes,&its);
86: SNESGetLinearSolveIterations(snes,&lits);
87: litspit = ((PetscReal)lits)/((PetscReal)its);
88: PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);
89: PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits);
90: PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / Newton = %e\n",litspit);
92: DMMGDestroy(dmmg);
93: PetscFinalize();
95: return 0;
96: }
97: /* -------------------- Evaluate Function F(x) --------------------- */
100: PetscErrorCode FormFunction(SNES snes,Vec T,Vec F,void* ptr)
101: {
102: DMMG dmmg = (DMMG)ptr;
104: PetscInt i,j,mx,my,xs,ys,xm,ym;
105: PetscScalar hx,hy;
106: PetscScalar **t,**f,gradup,graddown,gradleft,gradright,gradx,grady;
107: PetscScalar coeffup,coeffdown,coeffleft,coeffright;
108: Vec localT;
111: DAGetLocalVector((DA)dmmg->dm,&localT);
112: DAGetInfo((DA)dmmg->dm,PETSC_NULL,&mx,&my,0,0,0,0,0,0,0,0);
113: hx = 1.0/(PetscReal)(mx-1); hy = 1.0/(PetscReal)(my-1);
114:
115: /* Get ghost points */
116: DAGlobalToLocalBegin((DA)dmmg->dm,T,INSERT_VALUES,localT);
117: DAGlobalToLocalEnd((DA)dmmg->dm,T,INSERT_VALUES,localT);
118: DAGetCorners((DA)dmmg->dm,&xs,&ys,0,&xm,&ym,0);
119: DAVecGetArray((DA)dmmg->dm,localT,&t);
120: DAVecGetArray((DA)dmmg->dm,F,&f);
122: /* Evaluate function */
123: for (j=ys; j<ys+ym; j++) {
124: for (i=xs; i<xs+xm; i++) {
126: if (i == 0 || i == mx-1 || j == 0 || j == my-1) {
128: f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
129:
130: } else {
132: gradup = (t[j+1][i] - t[j][i])/hy;
133: graddown = (t[j][i] - t[j-1][i])/hy;
134: gradright = (t[j][i+1] - t[j][i])/hx;
135: gradleft = (t[j][i] - t[j][i-1])/hx;
137: gradx = .5*(t[j][i+1] - t[j][i-1])/hx;
138: grady = .5*(t[j+1][i] - t[j-1][i])/hy;
140: coeffup = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
141: coeffdown = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);
143: coeffleft = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
144: coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);
146: f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
147:
148: }
150: }
151: }
152: DAVecRestoreArray((DA)dmmg->dm,localT,&t);
153: DAVecRestoreArray((DA)dmmg->dm,F,&f);
154: DARestoreLocalVector((DA)dmmg->dm,&localT);
155: return(0);
156: }
158: PetscErrorCode FormFunctionLocal(DALocalInfo *info,PetscScalar **t,PetscScalar **f,void *ptr)
159: {
160: PetscInt i,j;
161: PetscScalar hx,hy;
162: PetscScalar gradup,graddown,gradleft,gradright,gradx,grady;
163: PetscScalar coeffup,coeffdown,coeffleft,coeffright;
166: hx = 1.0/(PetscReal)(info->mx-1); hy = 1.0/(PetscReal)(info->my-1);
167:
168: /* Evaluate function */
169: for (j=info->ys; j<info->ys+info->ym; j++) {
170: for (i=info->xs; i<info->xs+info->xm; i++) {
172: if (i == 0 || i == info->mx-1 || j == 0 || j == info->my-1) {
174: f[j][i] = t[j][i] - (1.0 - (2.0*hx*(PetscReal)i - 1.0)*(2.0*hx*(PetscReal)i - 1.0));
175:
176: } else {
178: gradup = (t[j+1][i] - t[j][i])/hy;
179: graddown = (t[j][i] - t[j-1][i])/hy;
180: gradright = (t[j][i+1] - t[j][i])/hx;
181: gradleft = (t[j][i] - t[j][i-1])/hx;
183: gradx = .5*(t[j][i+1] - t[j][i-1])/hx;
184: grady = .5*(t[j+1][i] - t[j-1][i])/hy;
186: coeffup = 1.0/PetscSqrtScalar(1.0 + gradup*gradup + gradx*gradx);
187: coeffdown = 1.0/PetscSqrtScalar(1.0 + graddown*graddown + gradx*gradx);
189: coeffleft = 1.0/PetscSqrtScalar(1.0 + gradleft*gradleft + grady*grady);
190: coeffright = 1.0/PetscSqrtScalar(1.0 + gradright*gradright + grady*grady);
192: f[j][i] = (coeffup*gradup - coeffdown*graddown)*hx + (coeffright*gradright - coeffleft*gradleft)*hy;
193:
194: }
196: }
197: }
198: return(0);
199: }