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: /*T
 10:    Concepts: SNES^sequential Bratu example
 11:    Processors: 1
 12: T*/

 14: /* ------------------------------------------------------------------------

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

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

 21:     with boundary conditions

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

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

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

 31:   ------------------------------------------------------------------------- */

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

 43: #include <petscsnes.h>

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

 56: /*
 57:    User-defined routines
 58: */
 59: extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 60: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 61: extern PetscErrorCode FormInitialGuess(AppCtx*,Vec);
 62: extern PetscErrorCode ConvergenceTest(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
 63: extern PetscErrorCode ConvergenceDestroy(void*);
 64: extern PetscErrorCode postcheck(SNES,Vec,Vec,Vec,PetscBool*,PetscBool*,void*);

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

 81:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 82:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 83:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");

 85:   /*
 86:      Initialize problem parameters
 87:   */
 88:   user.mx = 4; user.my = 4; user.param = 6.0;
 89:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
 90:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
 91:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
 92:   PetscOptionsGetBool(NULL,NULL,"-pc",&pc,NULL);
 93:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lambda is out of range");
 94:   N = user.mx*user.my;
 95:   PetscOptionsGetBool(NULL,NULL,"-use_convergence_test",&use_convergence_test,NULL);

 97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 98:      Create nonlinear solver context
 99:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

101:   SNESCreate(PETSC_COMM_WORLD,&snes);

103:   if (pc) {
104:     SNESSetType(snes,SNESNEWTONTR);
105:     SNESNewtonTRSetPostCheck(snes, postcheck,NULL);
106:   }

108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109:      Create vector data structures; set function evaluation routine
110:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

112:   VecCreate(PETSC_COMM_WORLD,&x);
113:   VecSetSizes(x,PETSC_DECIDE,N);
114:   VecSetFromOptions(x);
115:   VecDuplicate(x,&r);

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

127:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128:      Create matrix data structure; set Jacobian evaluation routine
129:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

146:   /*
147:      This option will cause the Jacobian to be computed via finite differences
148:     efficiently using a coloring of the columns of the matrix.
149:   */
150:   PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
151:   if (matrix_free && fd_coloring) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_INCOMP,"Use only one of -snes_mf, -snes_fd_coloring options!\nYou can do -snes_mf_operator -snes_fd_coloring");

153:   if (fd_coloring) {
154:     ISColoring   iscoloring;
155:     MatColoring  mc;

157:     /*
158:       This initializes the nonzero structure of the Jacobian. This is artificial
159:       because clearly if we had a routine to compute the Jacobian we won't need
160:       to use finite differences.
161:     */
162:     FormJacobian(snes,x,J,J,&user);

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

207:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208:      Customize nonlinear solver; set runtime options
209:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

211:   /*
212:      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
213:   */
214:   SNESSetFromOptions(snes);

216:   /*
217:      Set array that saves the function norms.  This array is intended
218:      when the user wants to save the convergence history for later use
219:      rather than just to view the function norms via -snes_monitor.
220:   */
221:   SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);

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

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

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

259:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260:      Free work space.  All PETSc objects should be destroyed when they
261:      are no longer needed.
262:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

264:   if (!matrix_free) {
265:     MatDestroy(&J);
266:   }
267:   if (fd_coloring) {
268:     MatFDColoringDestroy(&fdcoloring);
269:   }
270:   VecDestroy(&x);
271:   VecDestroy(&r);
272:   SNESDestroy(&snes);
273:   PetscFinalize();
274:   return ierr;
275: }
276: /* ------------------------------------------------------------------- */
277: /*
278:    FormInitialGuess - Forms initial approximation.

280:    Input Parameters:
281:    user - user-defined application context
282:    X - vector

284:    Output Parameter:
285:    X - vector
286:  */
287: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
288: {
289:   PetscInt       i,j,row,mx,my;
291:   PetscReal      lambda,temp1,temp,hx,hy;
292:   PetscScalar    *x;

294:   mx     = user->mx;
295:   my     = user->my;
296:   lambda = user->param;

298:   hx = 1.0 / (PetscReal)(mx-1);
299:   hy = 1.0 / (PetscReal)(my-1);

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

322:   /*
323:      Restore vector
324:   */
325:   VecRestoreArray(X,&x);
326:   return 0;
327: }
328: /* ------------------------------------------------------------------- */
329: /*
330:    FormFunction - Evaluates nonlinear function, F(x).

332:    Input Parameters:
333: .  snes - the SNES context
334: .  X - input vector
335: .  ptr - optional user-defined context, as set by SNESSetFunction()

337:    Output Parameter:
338: .  F - function vector
339:  */
340: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
341: {
342:   AppCtx            *user = (AppCtx*)ptr;
343:   PetscInt          i,j,row,mx,my;
344:   PetscErrorCode    ierr;
345:   PetscReal         two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx;
346:   PetscScalar       ut,ub,ul,ur,u,uxx,uyy,sc,*f;
347:   const PetscScalar *x;

349:   mx     = user->mx;
350:   my     = user->my;
351:   lambda = user->param;
352:   hx     = one / (PetscReal)(mx-1);
353:   hy     = one / (PetscReal)(my-1);
354:   sc     = hx*hy;
355:   hxdhy  = hx/hy;
356:   hydhx  = hy/hx;

358:   /*
359:      Get pointers to vector data
360:   */
361:   VecGetArrayRead(X,&x);
362:   VecGetArray(F,&f);

364:   /*
365:      Compute function
366:   */
367:   for (j=0; j<my; j++) {
368:     for (i=0; i<mx; i++) {
369:       row = i + j*mx;
370:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
371:         f[row] = x[row];
372:         continue;
373:       }
374:       u      = x[row];
375:       ub     = x[row - mx];
376:       ul     = x[row - 1];
377:       ut     = x[row + mx];
378:       ur     = x[row + 1];
379:       uxx    = (-ur + two*u - ul)*hydhx;
380:       uyy    = (-ut + two*u - ub)*hxdhy;
381:       f[row] = uxx + uyy - sc*lambda*PetscExpScalar(u);
382:     }
383:   }

385:   /*
386:      Restore vectors
387:   */
388:   VecRestoreArrayRead(X,&x);
389:   VecRestoreArray(F,&f);
390:   return 0;
391: }
392: /* ------------------------------------------------------------------- */
393: /*
394:    FormJacobian - Evaluates Jacobian matrix.

396:    Input Parameters:
397: .  snes - the SNES context
398: .  x - input vector
399: .  ptr - optional user-defined context, as set by SNESSetJacobian()

401:    Output Parameters:
402: .  A - Jacobian matrix
403: .  B - optionally different preconditioning matrix
404: .  flag - flag indicating matrix structure
405: */
406: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
407: {
408:   AppCtx            *user = (AppCtx*)ptr;   /* user-defined applicatin context */
409:   PetscInt          i,j,row,mx,my,col[5];
410:   PetscErrorCode    ierr;
411:   PetscScalar       two = 2.0,one = 1.0,lambda,v[5],sc;
412:   const PetscScalar *x;
413:   PetscReal         hx,hy,hxdhy,hydhx;

415:   mx     = user->mx;
416:   my     = user->my;
417:   lambda = user->param;
418:   hx     = 1.0 / (PetscReal)(mx-1);
419:   hy     = 1.0 / (PetscReal)(my-1);
420:   sc     = hx*hy;
421:   hxdhy  = hx/hy;
422:   hydhx  = hy/hx;

424:   /*
425:      Get pointer to vector data
426:   */
427:   VecGetArrayRead(X,&x);

429:   /*
430:      Compute entries of the Jacobian
431:   */
432:   for (j=0; j<my; j++) {
433:     for (i=0; i<mx; i++) {
434:       row = i + j*mx;
435:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
436:         MatSetValues(jac,1,&row,1,&row,&one,INSERT_VALUES);
437:         continue;
438:       }
439:       v[0] = -hxdhy; col[0] = row - mx;
440:       v[1] = -hydhx; col[1] = row - 1;
441:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = row;
442:       v[3] = -hydhx; col[3] = row + 1;
443:       v[4] = -hxdhy; col[4] = row + mx;
444:       MatSetValues(jac,1,&row,5,col,v,INSERT_VALUES);
445:     }
446:   }

448:   /*
449:      Restore vector
450:   */
451:   VecRestoreArrayRead(X,&x);

453:   /*
454:      Assemble matrix
455:   */
456:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
457:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

459:   if (jac != J) {
460:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
461:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
462:   }

464:   return 0;
465: }

467: PetscErrorCode ConvergenceTest(KSP ksp,PetscInt it,PetscReal nrm,KSPConvergedReason *reason,void *ctx)
468: {

472:   *reason = KSP_CONVERGED_ITERATING;
473:   if (it > 1) {
474:     *reason = KSP_CONVERGED_ITS;
475:     PetscInfo(NULL,"User provided convergence test returning after 2 iterations\n");
476:   }
477:   return(0);
478: }

480: PetscErrorCode ConvergenceDestroy(void* ctx)
481: {

485:   PetscInfo(NULL,"User provided convergence destroy called\n");
486:   PetscFree(ctx);
487:   return(0);
488: }

490: PetscErrorCode postcheck(SNES snes,Vec x,Vec y,Vec w,PetscBool *changed_y,PetscBool *changed_w,void *ctx)
491: {
493:   PetscReal      norm;
494:   Vec            tmp;

497:   VecDuplicate(x,&tmp);
498:   VecWAXPY(tmp,-1.0,x,w);
499:   VecNorm(tmp,NORM_2,&norm);
500:   VecDestroy(&tmp);
501:   PetscPrintf(PETSC_COMM_WORLD,"Norm of search step %g\n",(double)norm);
502:   return(0);
503: }

505: /*TEST

507:    build:
508:       requires: !single

510:    test:
511:       args: -ksp_gmres_cgs_refinement_type refine_always

513:    test:
514:       suffix: 2
515:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always

517:    test:
518:       suffix: 2a
519:       filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
520:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
521:       requires: defined(PETSC_USE_INFO)

523:    test:
524:       suffix: 2b
525:       filter: grep -i  "User provided convergence test" > /dev/null  && echo "Found User provided convergence test"
526:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -use_convergence_test -info
527:       requires: defined(PETSC_USE_INFO)

529:    test:
530:       suffix: 3
531:       args: -snes_monitor_short -mat_coloring_type sl -snes_fd_coloring -mx 8 -my 11 -ksp_gmres_cgs_refinement_type refine_always

533:    test:
534:       suffix: 4
535:       args: -pc -par 6.807 -snes_monitor -snes_converged_reason

537: TEST*/