Actual source code: ex1.c


  2: static char help[] = "Solves the nonlinear system, the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular domain.\n\
  3: This example also illustrates the use of matrix coloring.  Runtime options include:\n\
  4:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  5:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
  6:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
  7:   -my <yg>, where <yg> = number of grid points in the y-direction\n\n";

  9: /* ------------------------------------------------------------------------

 11:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 12:     the partial differential equation

 14:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,

 16:     with boundary conditions

 18:              u = 0  for  x = 0, x = 1, y = 0, y = 1.

 20:     A finite difference approximation with the usual 5-point stencil
 21:     is used to discretize the boundary value problem to obtain a nonlinear
 22:     system of equations.

 24:     The parallel version of this code is snes/tutorials/ex5.c

 26:   ------------------------------------------------------------------------- */

 28: /*
 29:    Include "petscsnes.h" so that we can use SNES solvers.  Note that
 30:    this file automatically includes:
 31:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 32:      petscmat.h - matrices
 33:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 34:      petscviewer.h - viewers               petscpc.h  - preconditioners
 35:      petscksp.h   - linear solvers
 36: */

 38: #include <petscsnes.h>

 40: /*
 41:    User-defined application context - contains data needed by the
 42:    application-provided call-back routines, FormJacobian() and
 43:    FormFunction().
 44: */
 45: typedef struct {
 46:   PetscReal param;              /* test problem parameter */
 47:   PetscInt  mx;                 /* Discretization in x-direction */
 48:   PetscInt  my;                 /* Discretization in y-direction */
 49: } AppCtx;

 51: /*
 52:    User-defined routines
 53: */
 54: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 55: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 56: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
 57: extern PetscErrorCode ConvergenceTest(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
 58: extern PetscErrorCode ConvergenceDestroy(void*);
 59: extern PetscErrorCode postcheck(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);

 61: int main(int argc,char **argv)
 62: {
 63:   SNES           snes;                 /* nonlinear solver context */
 64:   Vec            x,r;                 /* solution, residual vectors */
 65:   Mat            J;                    /* Jacobian matrix */
 66:   AppCtx         user;                 /* user-defined application context */
 67:   PetscInt       i,its,N,hist_its[50];
 68:   PetscMPIInt    size;
 69:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.,history[50];
 70:   MatFDColoring  fdcoloring;
 71:   PetscBool      matrix_free = PETSC_FALSE,flg,fd_coloring = PETSC_FALSE, use_convergence_test = PETSC_FALSE,pc = PETSC_FALSE;
 72:   KSP            ksp;
 73:   PetscInt       *testarray;

 75:   PetscInitialize(&argc,&argv,(char*)0,help);
 76:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 79:   /*
 80:      Initialize problem parameters
 81:   */
 82:   user.mx = 4; user.my = 4; user.param = 6.0;
 83:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
 84:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
 85:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
 86:   PetscOptionsGetBool(NULL,NULL,"-pc",&pc,NULL);
 88:   N = user.mx*user.my;
 89:   PetscOptionsGetBool(NULL,NULL,"-use_convergence_test",&use_convergence_test,NULL);

 91:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 92:      Create nonlinear solver context
 93:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 95:   SNESCreate(PETSC_COMM_WORLD,&snes);

 97:   if (pc) {
 98:     SNESSetType(snes,SNESNEWTONTR);
 99:     SNESNewtonTRSetPostCheck(snes, postcheck,NULL);
100:   }

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:      Create vector data structures; set function evaluation routine
104:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   VecCreate(PETSC_COMM_WORLD,&x);
107:   VecSetSizes(x,PETSC_DECIDE,N);
108:   VecSetFromOptions(x);
109:   VecDuplicate(x,&r);

111:   /*
112:      Set function evaluation routine and vector.  Whenever the nonlinear
113:      solver needs to evaluate the nonlinear function, it will call this
114:      routine.
115:       - Note that the final routine argument is the user-defined
116:         context that provides application-specific data for the
117:         function evaluation routine.
118:   */
119:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

121:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122:      Create matrix data structure; set Jacobian evaluation routine
123:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

125:   /*
126:      Create matrix. Here we only approximately preallocate storage space
127:      for the Jacobian.  See the users manual for a discussion of better
128:      techniques for preallocating matrix memory.
129:   */
130:   PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);
131:   if (!matrix_free) {
132:     PetscBool matrix_free_operator = PETSC_FALSE;
133:     PetscOptionsGetBool(NULL,NULL,"-snes_mf_operator",&matrix_free_operator,NULL);
134:     if (matrix_free_operator) matrix_free = PETSC_FALSE;
135:   }
136:   if (!matrix_free) {
137:     MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
138:   }

140:   /*
141:      This option will cause the Jacobian to be computed via finite differences
142:     efficiently using a coloring of the columns of the matrix.
143:   */
144:   PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);

147:   if (fd_coloring) {
148:     ISColoring   iscoloring;
149:     MatColoring  mc;

151:     /*
152:       This initializes the nonzero structure of the Jacobian. This is artificial
153:       because clearly if we had a routine to compute the Jacobian we won't need
154:       to use finite differences.
155:     */
156:     FormJacobian(snes,x,J,J,&user);

158:     /*
159:        Color the matrix, i.e. determine groups of columns that share no common
160:       rows. These columns in the Jacobian can all be computed simultaneously.
161:     */
162:     MatColoringCreate(J,&mc);
163:     MatColoringSetType(mc,MATCOLORINGSL);
164:     MatColoringSetFromOptions(mc);
165:     MatColoringApply(mc,&iscoloring);
166:     MatColoringDestroy(&mc);
167:     /*
168:        Create the data structure that SNESComputeJacobianDefaultColor() uses
169:        to compute the actual Jacobians via finite differences.
170:     */
171:     MatFDColoringCreate(J,iscoloring,&fdcoloring);
172:     MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
173:     MatFDColoringSetFromOptions(fdcoloring);
174:     MatFDColoringSetUp(J,iscoloring,fdcoloring);
175:     /*
176:         Tell SNES to use the routine SNESComputeJacobianDefaultColor()
177:       to compute Jacobians.
178:     */
179:     SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,fdcoloring);
180:     ISColoringDestroy(&iscoloring);
181:   }
182:   /*
183:      Set Jacobian matrix data structure and default Jacobian evaluation
184:      routine.  Whenever the nonlinear solver needs to compute the
185:      Jacobian matrix, it will call this routine.
186:       - Note that the final routine argument is the user-defined
187:         context that provides application-specific data for the
188:         Jacobian evaluation routine.
189:       - The user can override with:
190:          -snes_fd : default finite differencing approximation of Jacobian
191:          -snes_mf : matrix-free Newton-Krylov method with no preconditioning
192:                     (unless user explicitly sets preconditioner)
193:          -snes_mf_operator : form preconditioning matrix as set by the user,
194:                              but use matrix-free approx for Jacobian-vector
195:                              products within Newton-Krylov method
196:   */
197:   else if (!matrix_free) {
198:     SNESSetJacobian(snes,J,J,FormJacobian,(void*)&user);
199:   }

201:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202:      Customize nonlinear solver; set runtime options
203:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

205:   /*
206:      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
207:   */
208:   SNESSetFromOptions(snes);

210:   /*
211:      Set array that saves the function norms.  This array is intended
212:      when the user wants to save the convergence history for later use
213:      rather than just to view the function norms via -snes_monitor.
214:   */
215:   SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);

217:   /*
218:       Add a user provided convergence test; this is to test that SNESNEWTONTR properly calls the
219:       user provided test before the specialized test. The convergence context is just an array to
220:       test that it gets properly freed at the end
221:   */
222:   if (use_convergence_test) {
223:     SNESGetKSP(snes,&ksp);
224:     PetscMalloc1(5,&testarray);
225:     KSPSetConvergenceTest(ksp,ConvergenceTest,testarray,ConvergenceDestroy);
226:   }

228:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229:      Evaluate initial guess; then solve nonlinear system
230:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231:   /*
232:      Note: The user should initialize the vector, x, with the initial guess
233:      for the nonlinear solver prior to calling SNESSolve().  In particular,
234:      to employ an initial guess of zero, the user should explicitly set
235:      this vector to zero by calling VecSet().
236:   */
237:   FormInitialGuess(&user,x);
238:   SNESSolve(snes,NULL,x);
239:   SNESGetIterationNumber(snes,&its);
240:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its);

242:   /*
243:      Print the convergence history.  This is intended just to demonstrate
244:      use of the data attained via SNESSetConvergenceHistory().
245:   */
246:   PetscOptionsHasName(NULL,NULL,"-print_history",&flg);
247:   if (flg) {
248:     for (i=0; i<its+1; i++) {
249:       PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %g\n",i,hist_its[i],(double)history[i]);
250:     }
251:   }

253:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254:      Free work space.  All PETSc objects should be destroyed when they
255:      are no longer needed.
256:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

258:   if (!matrix_free) {
259:     MatDestroy(&J);
260:   }
261:   if (fd_coloring) {
262:     MatFDColoringDestroy(&fdcoloring);
263:   }
264:   VecDestroy(&x);
265:   VecDestroy(&r);
266:   SNESDestroy(&snes);
267:   PetscFinalize();
268:   return 0;
269: }
270: /* ------------------------------------------------------------------- */
271: /*
272:    FormInitialGuess - Forms initial approximation.

274:    Input Parameters:
275:    user - user-defined application context
276:    X - vector

278:    Output Parameter:
279:    X - vector
280:  */
281: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
282: {
283:   PetscInt       i,j,row,mx,my;
284:   PetscReal      lambda,temp1,temp,hx,hy;
285:   PetscScalar    *x;

287:   mx     = user->mx;
288:   my     = user->my;
289:   lambda = user->param;

291:   hx = 1.0 / (PetscReal)(mx-1);
292:   hy = 1.0 / (PetscReal)(my-1);

294:   /*
295:      Get a pointer to vector data.
296:        - For default PETSc vectors, VecGetArray() returns a pointer to
297:          the data array.  Otherwise, the routine is implementation dependent.
298:        - You MUST call VecRestoreArray() when you no longer need access to
299:          the array.
300:   */
301:   VecGetArray(X,&x);
302:   temp1 = lambda/(lambda + 1.0);
303:   for (j=0; j<my; j++) {
304:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
305:     for (i=0; i<mx; i++) {
306:       row = i + j*mx;
307:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
308:         x[row] = 0.0;
309:         continue;
310:       }
311:       x[row] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
312:     }
313:   }

315:   /*
316:      Restore vector
317:   */
318:   VecRestoreArray(X,&x);
319:   return 0;
320: }
321: /* ------------------------------------------------------------------- */
322: /*
323:    FormFunction - Evaluates nonlinear function, F(x).

325:    Input Parameters:
326: .  snes - the SNES context
327: .  X - input vector
328: .  ptr - optional user-defined context, as set by SNESSetFunction()

330:    Output Parameter:
331: .  F - function vector
332:  */
333: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
334: {
335:   AppCtx            *user = (AppCtx*)ptr;
336:   PetscInt          i,j,row,mx,my;
337:   PetscReal         two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
338:   PetscScalar       ut,ub,ul,ur,u,uxx,uyy,sc,*f;
339:   const PetscScalar *x;

341:   mx     = user->mx;
342:   my     = user->my;
343:   lambda = user->param;
344:   hx     = one / (PetscReal)(mx-1);
345:   hy     = one / (PetscReal)(my-1);
346:   sc     = hx*hy;
347:   hxdhy  = hx/hy;
348:   hydhx  = hy/hx;

350:   /*
351:      Get pointers to vector data
352:   */
353:   VecGetArrayRead(X,&x);
354:   VecGetArray(F,&f);

356:   /*
357:      Compute function
358:   */
359:   for (j=0; j<my; j++) {
360:     for (i=0; i<mx; i++) {
361:       row = i + j*mx;
362:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
363:         f[row] = x[row];
364:         continue;
365:       }
366:       u      = x[row];
367:       ub     = x[row - mx];
368:       ul     = x[row - 1];
369:       ut     = x[row + mx];
370:       ur     = x[row + 1];
371:       uxx    = (-ur + two*u - ul)*hydhx;
372:       uyy    = (-ut + two*u - ub)*hxdhy;
373:       f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
374:     }
375:   }

377:   /*
378:      Restore vectors
379:   */
380:   VecRestoreArrayRead(X,&x);
381:   VecRestoreArray(F,&f);
382:   return 0;
383: }
384: /* ------------------------------------------------------------------- */
385: /*
386:    FormJacobian - Evaluates Jacobian matrix.

388:    Input Parameters:
389: .  snes - the SNES context
390: .  x - input vector
391: .  ptr - optional user-defined context, as set by SNESSetJacobian()

393:    Output Parameters:
394: .  A - Jacobian matrix
395: .  B - optionally different preconditioning matrix
396: .  flag - flag indicating matrix structure
397: */
398: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
399: {
400:   AppCtx            *user = (AppCtx*)ptr;   /* user-defined applicatin context */
401:   PetscInt          i,j,row,mx,my,col[5];
402:   PetscScalar       two = 2.0,one = 1.0,lambda,v[5],sc;
403:   const PetscScalar *x;
404:   PetscReal         hx,hy,hxdhy,hydhx;

406:   mx     = user->mx;
407:   my     = user->my;
408:   lambda = user->param;
409:   hx     = 1.0 / (PetscReal)(mx-1);
410:   hy     = 1.0 / (PetscReal)(my-1);
411:   sc     = hx*hy;
412:   hxdhy  = hx/hy;
413:   hydhx  = hy/hx;

415:   /*
416:      Get pointer to vector data
417:   */
418:   VecGetArrayRead(X,&x);

420:   /*
421:      Compute entries of the Jacobian
422:   */
423:   for (j=0; j<my; j++) {
424:     for (i=0; i<mx; i++) {
425:       row = i + j*mx;
426:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
427:         MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
428:         continue;
429:       }
430:       v[0] = -hxdhy; col[0] = row - mx;
431:       v[1] = -hydhx; col[1] = row - 1;
432:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
433:       v[3] = -hydhx; col[3] = row + 1;
434:       v[4] = -hxdhy; col[4] = row + mx;
435:       MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
436:     }
437:   }

439:   /*
440:      Restore vector
441:   */
442:   VecRestoreArrayRead(X,&x);

444:   /*
445:      Assemble matrix
446:   */
447:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
448:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

450:   if (jac != J) {
451:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
452:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
453:   }

455:   return 0;
456: }

458: PetscErrorCode ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason *reason,void *ctx)
459: {
460:   *reason = KSP_CONVERGED_ITERATING;
461:   if (it > 1) {
462:     *reason = KSP_CONVERGED_ITS;
463:     PetscInfo(NULL,"User provided convergence test returning after 2 iterations\n");
464:   }
465:   return 0;
466: }

468: PetscErrorCode ConvergenceDestroy(void* ctx)
469: {
470:   PetscInfo(NULL,"User provided convergence destroy called\n");
471:   PetscFree(ctx);
472:   return 0;
473: }

475: PetscErrorCode postcheck(SNES snes,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
476: {
477:   PetscReal      norm;
478:   Vec            tmp;

480:   VecDuplicate(x,&tmp);
481:   VecWAXPY(tmp,-1.0,x,w);
482:   VecNorm(tmp,NORM_2,&norm);
483:   VecDestroy(&tmp);
484:   PetscPrintf(PETSC_COMM_WORLD,"Norm of search step %g\n",(double)norm);
485:   return 0;
486: }

488: /*TEST

490:    build:
491:       requires: !single

493:    test:
494:       args: -ksp_gmres_cgs_refinement_type refine_always

496:    test:
497:       suffix: 2
498:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always

500:    test:
501:       suffix: 2a
502:       filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
503:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
504:       requires: defined(PETSC_USE_INFO)

506:    test:
507:       suffix: 2b
508:       filter: grep -i  "User provided convergence test" > /dev/null  && echo "Found User provided convergence test"
509:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
510:       requires: defined(PETSC_USE_INFO)

512:    test:
513:       suffix: 3
514:       args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always

516:    test:
517:       suffix: 4
518:       args: -pc -par 6.807 -snes_monitor -snes_converged_reason

520: TEST*/