Actual source code: ex35.c
1: static const char help[] = "-Laplacian u = b as a nonlinear problem.\n\n";
3: /*
5: The linear and nonlinear versions of these should give almost identical results on this problem
7: Richardson
8: Nonlinear:
9: -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor
11: Linear:
12: -snes_rtol 1.e-12 -snes_monitor -ksp_rtol 1.e-12 -ksp_monitor -ksp_type richardson -pc_type none -ksp_richardson_self_scale -info
14: GMRES
15: Nonlinear:
16: -snes_rtol 1.e-12 -snes_monitor -snes_type ngmres
18: Linear:
19: -snes_rtol 1.e-12 -snes_monitor -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
21: CG
22: Nonlinear:
23: -snes_rtol 1.e-12 -snes_monitor -snes_type ncg -snes_linesearch_monitor
25: Linear:
26: -snes_rtol 1.e-12 -snes_monitor -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none
28: Multigrid
29: Linear:
30: 1 level:
31: -snes_rtol 1.e-12 -snes_monitor -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor
32: -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12 -ksp_monitor_true_residual
34: n levels:
35: -da_refine n
37: Nonlinear:
38: 1 level:
39: -snes_rtol 1.e-12 -snes_monitor -snes_type fas -fas_levels_snes_monitor
41: n levels:
42: -da_refine n -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly
44: */
46: /*
47: Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
48: Include "petscsnes.h" so that we can use SNES solvers. Note that this
49: */
50: #include <petscdm.h>
51: #include <petscdmda.h>
52: #include <petscsnes.h>
54: /*
55: User-defined routines
56: */
57: extern PetscErrorCode FormMatrix(DM,Mat);
58: extern PetscErrorCode MyComputeFunction(SNES,Vec,Vec,void*);
59: extern PetscErrorCode MyComputeJacobian(SNES,Vec,Mat,Mat,void*);
60: extern PetscErrorCode NonlinearGS(SNES,Vec);
62: int main(int argc,char **argv)
63: {
64: SNES snes; /* nonlinear solver */
65: SNES psnes; /* nonlinear Gauss-Seidel approximate solver */
66: Vec x,b; /* solution vector */
67: PetscInt its; /* iterations for convergence */
68: DM da;
69: PetscBool use_ngs_as_npc = PETSC_FALSE; /* use the nonlinear Gauss-Seidel approximate solver */
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Initialize program
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: PetscInitialize(&argc,&argv,(char*)0,help);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create nonlinear solver context
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: SNESCreate(PETSC_COMM_WORLD,&snes);
82: PetscOptionsGetBool(NULL,NULL,"-use_ngs_as_npc",&use_ngs_as_npc,0);
84: if (use_ngs_as_npc) {
85: SNESGetNPC(snes,&psnes);
86: SNESSetType(psnes,SNESSHELL);
87: SNESShellSetSolve(psnes,NonlinearGS);
88: }
90: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
91: Create distributed array (DMDA) to manage parallel grid and vectors
92: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
93: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
94: DMSetFromOptions(da);
95: DMSetUp(da);
96: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
97: SNESSetDM(snes,da);
98: if (use_ngs_as_npc) {
99: SNESShellSetContext(psnes,da);
100: }
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Extract global vectors from DMDA; then duplicate for remaining
103: vectors that are the same types
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: DMCreateGlobalVector(da,&x);
106: DMCreateGlobalVector(da,&b);
107: VecSet(b,1.0);
109: SNESSetFunction(snes,NULL,MyComputeFunction,NULL);
110: SNESSetJacobian(snes,NULL,NULL,MyComputeJacobian,NULL);
112: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113: Customize nonlinear solver; set runtime options
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115: SNESSetFromOptions(snes);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Solve nonlinear system
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: SNESSolve(snes,b,x);
121: SNESGetIterationNumber(snes,&its);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Free work space. All PETSc objects should be destroyed when they
128: are no longer needed.
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: VecDestroy(&x);
131: VecDestroy(&b);
132: SNESDestroy(&snes);
133: DMDestroy(&da);
134: PetscFinalize();
135: return 0;
136: }
138: /* ------------------------------------------------------------------- */
139: PetscErrorCode MyComputeFunction(SNES snes,Vec x,Vec F,void *ctx)
140: {
141: Mat J;
142: DM dm;
145: SNESGetDM(snes,&dm);
146: DMGetApplicationContext(dm,&J);
147: if (!J) {
148: DMSetMatType(dm,MATAIJ);
149: DMCreateMatrix(dm,&J);
150: MatSetDM(J, NULL);
151: FormMatrix(dm,J);
152: DMSetApplicationContext(dm,J);
153: DMSetApplicationContextDestroy(dm,(PetscErrorCode (*)(void**))MatDestroy);
154: }
155: MatMult(J,x,F);
156: return 0;
157: }
159: PetscErrorCode MyComputeJacobian(SNES snes,Vec x,Mat J,Mat Jp,void *ctx)
160: {
161: DM dm;
164: SNESGetDM(snes,&dm);
165: FormMatrix(dm,Jp);
166: return 0;
167: }
169: PetscErrorCode FormMatrix(DM da,Mat jac)
170: {
171: PetscInt i,j,nrows = 0;
172: MatStencil col[5],row,*rows;
173: PetscScalar v[5],hx,hy,hxdhy,hydhx;
174: DMDALocalInfo info;
177: DMDAGetLocalInfo(da,&info);
178: hx = 1.0/(PetscReal)(info.mx-1);
179: hy = 1.0/(PetscReal)(info.my-1);
180: hxdhy = hx/hy;
181: hydhx = hy/hx;
183: PetscMalloc1(info.ym*info.xm,&rows);
184: /*
185: Compute entries for the locally owned part of the Jacobian.
186: - Currently, all PETSc parallel matrix formats are partitioned by
187: contiguous chunks of rows across the processors.
188: - Each processor needs to insert only elements that it owns
189: locally (but any non-local elements will be sent to the
190: appropriate processor during matrix assembly).
191: - Here, we set all entries for a particular row at once.
192: - We can set matrix entries either using either
193: MatSetValuesLocal() or MatSetValues(), as discussed above.
194: */
195: for (j=info.ys; j<info.ys+info.ym; j++) {
196: for (i=info.xs; i<info.xs+info.xm; i++) {
197: row.j = j; row.i = i;
198: /* boundary points */
199: if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
200: v[0] = 2.0*(hydhx + hxdhy);
201: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
202: rows[nrows].i = i;
203: rows[nrows++].j = j;
204: } else {
205: /* interior grid points */
206: v[0] = -hxdhy; col[0].j = j - 1; col[0].i = i;
207: v[1] = -hydhx; col[1].j = j; col[1].i = i-1;
208: v[2] = 2.0*(hydhx + hxdhy); col[2].j = row.j; col[2].i = row.i;
209: v[3] = -hydhx; col[3].j = j; col[3].i = i+1;
210: v[4] = -hxdhy; col[4].j = j + 1; col[4].i = i;
211: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
212: }
213: }
214: }
216: /*
217: Assemble matrix, using the 2-step process:
218: MatAssemblyBegin(), MatAssemblyEnd().
219: */
220: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
221: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
222: MatZeroRowsColumnsStencil(jac,nrows,rows,2.0*(hydhx + hxdhy),NULL,NULL);
223: PetscFree(rows);
224: /*
225: Tell the matrix we will never add a new nonzero location to the
226: matrix. If we do, it will generate an error.
227: */
228: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
229: return 0;
230: }
232: /* ------------------------------------------------------------------- */
233: /*
234: Applies some sweeps on nonlinear Gauss-Seidel on each process
236: */
237: PetscErrorCode NonlinearGS(SNES snes,Vec X)
238: {
239: PetscInt i,j,Mx,My,xs,ys,xm,ym,its,l;
240: PetscReal hx,hy,hxdhy,hydhx;
241: PetscScalar **x,F,J,u,uxx,uyy;
242: DM da;
243: Vec localX;
246: SNESGetTolerances(snes,NULL,NULL,NULL,&its,NULL);
247: SNESShellGetContext(snes,&da);
249: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
251: hx = 1.0/(PetscReal)(Mx-1);
252: hy = 1.0/(PetscReal)(My-1);
253: hxdhy = hx/hy;
254: hydhx = hy/hx;
256: DMGetLocalVector(da,&localX);
258: for (l=0; l<its; l++) {
260: DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
261: DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
262: /*
263: Get a pointer to vector data.
264: - For default PETSc vectors, VecGetArray() returns a pointer to
265: the data array. Otherwise, the routine is implementation dependent.
266: - You MUST call VecRestoreArray() when you no longer need access to
267: the array.
268: */
269: DMDAVecGetArray(da,localX,&x);
271: /*
272: Get local grid boundaries (for 2-dimensional DMDA):
273: xs, ys - starting grid indices (no ghost points)
274: xm, ym - widths of local grid (no ghost points)
276: */
277: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
279: for (j=ys; j<ys+ym; j++) {
280: for (i=xs; i<xs+xm; i++) {
281: if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
282: /* boundary conditions are all zero Dirichlet */
283: x[j][i] = 0.0;
284: } else {
285: u = x[j][i];
286: uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
287: uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
288: F = uxx + uyy;
289: J = 2.0*(hydhx + hxdhy);
290: u = u - F/J;
292: x[j][i] = u;
293: }
294: }
295: }
297: /*
298: Restore vector
299: */
300: DMDAVecRestoreArray(da,localX,&x);
301: DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
302: DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
303: }
304: DMRestoreLocalVector(da,&localX);
305: return 0;
306: }
308: /*TEST
310: test:
311: args: -snes_monitor_short -snes_type nrichardson
312: requires: !single
314: test:
315: suffix: 2
316: args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale
317: requires: !single
319: test:
320: suffix: 3
321: args: -snes_monitor_short -snes_type ngmres
323: test:
324: suffix: 4
325: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none
327: test:
328: suffix: 5
329: args: -snes_monitor_short -snes_type ncg
331: test:
332: suffix: 6
333: args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none
335: test:
336: suffix: 7
337: args: -da_refine 2 -snes_monitor_short -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor_short -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor_short
338: requires: !single
340: test:
341: suffix: 8
342: args: -da_refine 2 -snes_monitor_short -snes_type fas -fas_levels_snes_monitor_short -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly -snes_type fas -snes_rtol 1.e-5
344: test:
345: suffix: 9
346: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc
348: test:
349: suffix: 10
350: args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none -snes_type newtontrdc -snes_trdc_use_cauchy false
352: TEST*/