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;

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

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

 94:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 95:      Create nonlinear solver context
 96:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 98:   SNESCreate(PETSC_COMM_WORLD, &snes);

100:   if (pc) {
101:     SNESSetType(snes, SNESNEWTONTR);
102:     SNESNewtonTRSetPostCheck(snes, postcheck, NULL);
103:   }

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Create vector data structures; set function evaluation routine
107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

109:   VecCreate(PETSC_COMM_WORLD, &x);
110:   VecSetSizes(x, PETSC_DECIDE, N);
111:   VecSetFromOptions(x);
112:   VecDuplicate(x, &r);

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

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:      Create matrix data structure; set Jacobian evaluation routine
126:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

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

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

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

257:   if (!matrix_free) MatDestroy(&J);
258:   if (fd_coloring) MatFDColoringDestroy(&fdcoloring);
259:   VecDestroy(&x);
260:   VecDestroy(&r);
261:   SNESDestroy(&snes);
262:   PetscFinalize();
263:   return 0;
264: }

266: /*
267:    FormInitialGuess - Forms initial approximation.

269:    Input Parameters:
270:    user - user-defined application context
271:    X - vector

273:    Output Parameter:
274:    X - vector
275:  */
276: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
277: {
278:   PetscInt     i, j, row, mx, my;
279:   PetscReal    lambda, temp1, temp, hx, hy;
280:   PetscScalar *x;

282:   mx     = user->mx;
283:   my     = user->my;
284:   lambda = user->param;

286:   hx = 1.0 / (PetscReal)(mx - 1);
287:   hy = 1.0 / (PetscReal)(my - 1);

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

310:   /*
311:      Restore vector
312:   */
313:   VecRestoreArray(X, &x);
314:   return 0;
315: }

317: /*
318:    FormFunction - Evaluates nonlinear function, F(x).

320:    Input Parameters:
321: .  snes - the SNES context
322: .  X - input vector
323: .  ptr - optional user-defined context, as set by SNESSetFunction()

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

336:   mx     = user->mx;
337:   my     = user->my;
338:   lambda = user->param;
339:   hx     = one / (PetscReal)(mx - 1);
340:   hy     = one / (PetscReal)(my - 1);
341:   sc     = hx * hy;
342:   hxdhy  = hx / hy;
343:   hydhx  = hy / hx;

345:   /*
346:      Get pointers to vector data
347:   */
348:   VecGetArrayRead(X, &x);
349:   VecGetArray(F, &f);

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

372:   /*
373:      Restore vectors
374:   */
375:   VecRestoreArrayRead(X, &x);
376:   VecRestoreArray(F, &f);
377:   return 0;
378: }

380: /*
381:    FormJacobian - Evaluates Jacobian matrix.

383:    Input Parameters:
384: .  snes - the SNES context
385: .  x - input vector
386: .  ptr - optional user-defined context, as set by SNESSetJacobian()

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

401:   mx     = user->mx;
402:   my     = user->my;
403:   lambda = user->param;
404:   hx     = 1.0 / (PetscReal)(mx - 1);
405:   hy     = 1.0 / (PetscReal)(my - 1);
406:   sc     = hx * hy;
407:   hxdhy  = hx / hy;
408:   hydhx  = hy / hx;

410:   /*
411:      Get pointer to vector data
412:   */
413:   VecGetArrayRead(X, &x);

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