Actual source code: ex35.c

  1: static const char help[] = "-Laplacian u = b as a nonlinear problem.\n\n";

  3: /*T
  4:    Concepts: SNES^parallel Bratu example
  5:    Concepts: DMDA^using distributed arrays;
  6:    Concepts: IS coloirng types;
  7:    Processors: n
  8: T*/



 12: /*

 14:     The linear and nonlinear versions of these should give almost identical results on this problem

 16:     Richardson
 17:       Nonlinear:
 18:         -snes_rtol 1.e-12 -snes_monitor -snes_type nrichardson -snes_linesearch_monitor

 20:       Linear:
 21:         -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

 23:     GMRES
 24:       Nonlinear:
 25:        -snes_rtol 1.e-12 -snes_monitor  -snes_type ngmres

 27:       Linear:
 28:        -snes_rtol 1.e-12 -snes_monitor  -ksp_type gmres -ksp_monitor -ksp_rtol 1.e-12 -pc_type none

 30:     CG
 31:        Nonlinear:
 32:             -snes_rtol 1.e-12 -snes_monitor  -snes_type ncg -snes_linesearch_monitor

 34:        Linear:
 35:              -snes_rtol 1.e-12 -snes_monitor  -ksp_type cg -ksp_monitor -ksp_rtol 1.e-12 -pc_type none

 37:     Multigrid
 38:        Linear:
 39:           1 level:
 40:             -snes_rtol 1.e-12 -snes_monitor  -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type none -mg_levels_ksp_monitor
 41:             -mg_levels_ksp_richardson_self_scale -ksp_type richardson -ksp_monitor -ksp_rtol 1.e-12  -ksp_monitor_true_residual

 43:           n levels:
 44:             -da_refine n

 46:        Nonlinear:
 47:          1 level:
 48:            -snes_rtol 1.e-12 -snes_monitor  -snes_type fas -fas_levels_snes_monitor

 50:           n levels:
 51:             -da_refine n  -fas_coarse_snes_type newtonls -fas_coarse_pc_type lu -fas_coarse_ksp_type preonly

 53: */

 55: /*
 56:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
 57:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 58: */
 59: #include <petscdm.h>
 60: #include <petscdmda.h>
 61: #include <petscsnes.h>

 63: /*
 64:    User-defined routines
 65: */
 66: extern PetscErrorCode FormMatrix(DM,Mat);
 67: extern PetscErrorCode MyComputeFunction(SNES,Vec,Vec,void*);
 68: extern PetscErrorCode MyComputeJacobian(SNES,Vec,Mat,Mat,void*);
 69: extern PetscErrorCode NonlinearGS(SNES,Vec);

 71: int main(int argc,char **argv)
 72: {
 73:   SNES           snes;                                 /* nonlinear solver */
 74:   SNES           psnes;                                /* nonlinear Gauss-Seidel approximate solver */
 75:   Vec            x,b;                                  /* solution vector */
 76:   PetscInt       its;                                  /* iterations for convergence */
 78:   DM             da;
 79:   PetscBool      use_ngs_as_npc = PETSC_FALSE;                /* use the nonlinear Gauss-Seidel approximate solver */

 81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82:      Initialize program
 83:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 85:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88:      Create nonlinear solver context
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 90:   SNESCreate(PETSC_COMM_WORLD,&snes);

 92:   PetscOptionsGetBool(NULL,NULL,"-use_ngs_as_npc",&use_ngs_as_npc,0);

 94:   if (use_ngs_as_npc) {
 95:     SNESGetNPC(snes,&psnes);
 96:     SNESSetType(psnes,SNESSHELL);
 97:     SNESShellSetSolve(psnes,NonlinearGS);
 98:   }

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Create distributed array (DMDA) to manage parallel grid and vectors
102:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
104:   DMSetFromOptions(da);
105:   DMSetUp(da);
106:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
107:   SNESSetDM(snes,da);
108:   if (use_ngs_as_npc) {
109:     SNESShellSetContext(psnes,da);
110:   }
111:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112:      Extract global vectors from DMDA; then duplicate for remaining
113:      vectors that are the same types
114:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115:   DMCreateGlobalVector(da,&x);
116:   DMCreateGlobalVector(da,&b);
117:   VecSet(b,1.0);

119:   SNESSetFunction(snes,NULL,MyComputeFunction,NULL);
120:   SNESSetJacobian(snes,NULL,NULL,MyComputeJacobian,NULL);

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123:      Customize nonlinear solver; set runtime options
124:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125:   SNESSetFromOptions(snes);

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Solve nonlinear system
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130:   SNESSolve(snes,b,x);
131:   SNESGetIterationNumber(snes,&its);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

136:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137:      Free work space.  All PETSc objects should be destroyed when they
138:      are no longer needed.
139:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140:   VecDestroy(&x);
141:   VecDestroy(&b);
142:   SNESDestroy(&snes);
143:   DMDestroy(&da);
144:   PetscFinalize();
145:   return ierr;
146: }

148: /* ------------------------------------------------------------------- */
149: PetscErrorCode MyComputeFunction(SNES snes,Vec x,Vec F,void *ctx)
150: {
152:   Mat            J;
153:   DM             dm;

156:   SNESGetDM(snes,&dm);
157:   DMGetApplicationContext(dm,&J);
158:   if (!J) {
159:     DMSetMatType(dm,MATAIJ);
160:     DMCreateMatrix(dm,&J);
161:     MatSetDM(J, NULL);
162:     FormMatrix(dm,J);
163:     DMSetApplicationContext(dm,J);
164:     DMSetApplicationContextDestroy(dm,(PetscErrorCode (*)(void**))MatDestroy);
165:   }
166:   MatMult(J,x,F);
167:   return(0);
168: }

170: PetscErrorCode MyComputeJacobian(SNES snes,Vec x,Mat J,Mat Jp,void *ctx)
171: {
173:   DM             dm;

176:   SNESGetDM(snes,&dm);
177:   FormMatrix(dm,Jp);
178:   return(0);
179: }

181: PetscErrorCode FormMatrix(DM da,Mat jac)
182: {
184:   PetscInt       i,j,nrows = 0;
185:   MatStencil     col[5],row,*rows;
186:   PetscScalar    v[5],hx,hy,hxdhy,hydhx;
187:   DMDALocalInfo  info;

190:   DMDAGetLocalInfo(da,&info);
191:   hx    = 1.0/(PetscReal)(info.mx-1);
192:   hy    = 1.0/(PetscReal)(info.my-1);
193:   hxdhy = hx/hy;
194:   hydhx = hy/hx;

196:   PetscMalloc1(info.ym*info.xm,&rows);
197:   /*
198:      Compute entries for the locally owned part of the Jacobian.
199:       - Currently, all PETSc parallel matrix formats are partitioned by
200:         contiguous chunks of rows across the processors.
201:       - Each processor needs to insert only elements that it owns
202:         locally (but any non-local elements will be sent to the
203:         appropriate processor during matrix assembly).
204:       - Here, we set all entries for a particular row at once.
205:       - We can set matrix entries either using either
206:         MatSetValuesLocal() or MatSetValues(), as discussed above.
207:   */
208:   for (j=info.ys; j<info.ys+info.ym; j++) {
209:     for (i=info.xs; i<info.xs+info.xm; i++) {
210:       row.j = j; row.i = i;
211:       /* boundary points */
212:       if (i == 0 || j == 0 || i == info.mx-1 || j == info.my-1) {
213:         v[0]            = 2.0*(hydhx + hxdhy);
214:         MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
215:         rows[nrows].i   = i;
216:         rows[nrows++].j = j;
217:       } else {
218:         /* interior grid points */
219:         v[0] = -hxdhy;                                           col[0].j = j - 1; col[0].i = i;
220:         v[1] = -hydhx;                                           col[1].j = j;     col[1].i = i-1;
221:         v[2] = 2.0*(hydhx + hxdhy);                              col[2].j = row.j; col[2].i = row.i;
222:         v[3] = -hydhx;                                           col[3].j = j;     col[3].i = i+1;
223:         v[4] = -hxdhy;                                           col[4].j = j + 1; col[4].i = i;
224:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
225:       }
226:     }
227:   }

229:   /*
230:      Assemble matrix, using the 2-step process:
231:        MatAssemblyBegin(), MatAssemblyEnd().
232:   */
233:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
234:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
235:   MatZeroRowsColumnsStencil(jac,nrows,rows,2.0*(hydhx + hxdhy),NULL,NULL);
236:   PetscFree(rows);
237:   /*
238:      Tell the matrix we will never add a new nonzero location to the
239:      matrix. If we do, it will generate an error.
240:   */
241:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
242:   return(0);
243: }



247: /* ------------------------------------------------------------------- */
248: /*
249:       Applies some sweeps on nonlinear Gauss-Seidel on each process

251:  */
252: PetscErrorCode NonlinearGS(SNES snes,Vec X)
253: {
254:   PetscInt       i,j,Mx,My,xs,ys,xm,ym,its,l;
256:   PetscReal      hx,hy,hxdhy,hydhx;
257:   PetscScalar    **x,F,J,u,uxx,uyy;
258:   DM             da;
259:   Vec            localX;

262:   SNESGetTolerances(snes,NULL,NULL,NULL,&its,NULL);
263:   SNESShellGetContext(snes,(void**)&da);

265:   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);

267:   hx    = 1.0/(PetscReal)(Mx-1);
268:   hy    = 1.0/(PetscReal)(My-1);
269:   hxdhy = hx/hy;
270:   hydhx = hy/hx;


273:   DMGetLocalVector(da,&localX);

275:   for (l=0; l<its; l++) {

277:     DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
278:     DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);
279:     /*
280:      Get a pointer to vector data.
281:      - For default PETSc vectors, VecGetArray() returns a pointer to
282:      the data array.  Otherwise, the routine is implementation dependent.
283:      - You MUST call VecRestoreArray() when you no longer need access to
284:      the array.
285:      */
286:     DMDAVecGetArray(da,localX,&x);

288:     /*
289:      Get local grid boundaries (for 2-dimensional DMDA):
290:      xs, ys   - starting grid indices (no ghost points)
291:      xm, ym   - widths of local grid (no ghost points)

293:      */
294:     DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

296:     for (j=ys; j<ys+ym; j++) {
297:       for (i=xs; i<xs+xm; i++) {
298:         if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
299:           /* boundary conditions are all zero Dirichlet */
300:           x[j][i] = 0.0;
301:         } else {
302:           u   = x[j][i];
303:           uxx = (2.0*u - x[j][i-1] - x[j][i+1])*hydhx;
304:           uyy = (2.0*u - x[j-1][i] - x[j+1][i])*hxdhy;
305:           F   = uxx + uyy;
306:           J   = 2.0*(hydhx + hxdhy);
307:           u   = u - F/J;

309:           x[j][i] = u;
310:         }
311:       }
312:     }

314:     /*
315:      Restore vector
316:      */
317:     DMDAVecRestoreArray(da,localX,&x);
318:     DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X);
319:     DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X);
320:   }
321:   DMRestoreLocalVector(da,&localX);
322:   return(0);
323: }


326: /*TEST

328:    test:
329:       args: -snes_monitor_short -snes_type nrichardson
330:       requires: !single

332:    test:
333:       suffix: 2
334:       args: -snes_monitor_short -ksp_monitor_short -ksp_type richardson -pc_type none -ksp_richardson_self_scale
335:       requires: !single

337:    test:
338:       suffix: 3
339:       args: -snes_monitor_short -snes_type ngmres

341:    test:
342:       suffix: 4
343:       args: -snes_monitor_short -ksp_type gmres -ksp_monitor_short -pc_type none

345:    test:
346:       suffix: 5
347:       args: -snes_monitor_short -snes_type ncg

349:    test:
350:       suffix: 6
351:       args: -snes_monitor_short -ksp_type cg -ksp_monitor_short -pc_type none

353:    test:
354:       suffix: 7
355:       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
356:       requires: !single

358:    test:
359:       suffix: 8
360:       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

362: TEST*/