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*/