Actual source code: ex1.c

  2: /* Program usage:  ex4 [-help] [all PETSc options] */

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

 12: /*T
 13:    Concepts: SNES^sequential Bratu example
 14:    Processors: 1
 15: T*/

 17: /* ------------------------------------------------------------------------

 19:     Solid Fuel Ignition (SFI) problem.  This problem is modeled by
 20:     the partial differential equation
 21:   
 22:             -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
 23:   
 24:     with boundary conditions
 25:    
 26:              u = 0  for  x = 0, x = 1, y = 0, y = 1.
 27:   
 28:     A finite difference approximation with the usual 5-point stencil
 29:     is used to discretize the boundary value problem to obtain a nonlinear 
 30:     system of equations.

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

 34:   ------------------------------------------------------------------------- */

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

 46: #include <petscsnes.h>

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

 59: /* 
 60:    User-defined routines
 61: */

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

 81:   PetscInitialize(&argc,&argv,(char *)0,help);
 82:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 83:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"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(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
 90:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
 91:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
 92:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");
 93:   N = user.mx*user.my;
 94: 
 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:      Create nonlinear solver context
 97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 99:   SNESCreate(PETSC_COMM_WORLD,&snes);

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

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

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

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

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

139:   /*
140:      This option will cause the Jacobian to be computed via finite differences
141:     efficiently using a coloring of the columns of the matrix.
142:   */
143:   PetscOptionsGetBool(PETSC_NULL,"-snes_fd_coloring",&fd_coloring,PETSC_NULL);
144:   if (matrix_free && fd_coloring)  SETERRQ(PETSC_COMM_SELF,1,"Use only one of -snes_mf, -snes_fd_coloring options!\nYou can do -snes_mf_operator -snes_fd_coloring");

146:   if (fd_coloring) {
147:     ISColoring   iscoloring;
148:     MatStructure str;

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

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

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:      Customize nonlinear solver; set runtime options
197:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

199:   /*
200:      Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>)
201:   */
202:   SNESSetFromOptions(snes);

204:   /*
205:      Set array that saves the function norms.  This array is intended
206:      when the user wants to save the convergence history for later use
207:      rather than just to view the function norms via -snes_monitor.
208:   */
209:   SNESSetConvergenceHistory(snes,history,hist_its,50,PETSC_TRUE);

211:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
212:      Evaluate initial guess; then solve nonlinear system
213:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
214:   /*
215:      Note: The user should initialize the vector, x, with the initial guess
216:      for the nonlinear solver prior to calling SNESSolve().  In particular,
217:      to employ an initial guess of zero, the user should explicitly set
218:      this vector to zero by calling VecSet().
219:   */
220:   FormInitialGuess(&user,x);
221:   SNESSolve(snes,PETSC_NULL,x);
222:   SNESGetIterationNumber(snes,&its);
223:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);


226:   /* 
227:      Print the convergence history.  This is intended just to demonstrate
228:      use of the data attained via SNESSetConvergenceHistory().  
229:   */
230:   PetscOptionsHasName(PETSC_NULL,"-print_history",&flg);
231:   if (flg) {
232:     for (i=0; i<its+1; i++) {
233:       PetscPrintf(PETSC_COMM_WORLD,"iteration %D: Linear iterations %D Function norm = %G\n",i,hist_its[i],history[i]);
234:     }
235:   }

237:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
238:      Free work space.  All PETSc objects should be destroyed when they
239:      are no longer needed.
240:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

242:   if (!matrix_free) {
243:     MatDestroy(&J);
244:   }
245:   if (fd_coloring) {
246:     MatFDColoringDestroy(&fdcoloring);
247:   }
248:   VecDestroy(&x);
249:   VecDestroy(&r);
250:   SNESDestroy(&snes);
251:   PetscFinalize();

253:   return 0;
254: }
255: /* ------------------------------------------------------------------- */
258: /* 
259:    FormInitialGuess - Forms initial approximation.

261:    Input Parameters:
262:    user - user-defined application context
263:    X - vector

265:    Output Parameter:
266:    X - vector
267:  */
268: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
269: {
270:   PetscInt       i,j,row,mx,my;
272:   PetscReal      lambda,temp1,temp,hx,hy;
273:   PetscScalar    *x;

275:   mx         = user->mx;
276:   my         = user->my;
277:   lambda = user->param;

279:   hx    = 1.0 / (PetscReal)(mx-1);
280:   hy    = 1.0 / (PetscReal)(my-1);

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

303:   /*
304:      Restore vector
305:   */
306:   VecRestoreArray(X,&x);
307:   return 0;
308: }
309: /* ------------------------------------------------------------------- */
312: /* 
313:    FormFunction - Evaluates nonlinear function, F(x).

315:    Input Parameters:
316: .  snes - the SNES context
317: .  X - input vector
318: .  ptr - optional user-defined context, as set by SNESSetFunction()

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

331:   mx         = user->mx;
332:   my         = user->my;
333:   lambda = user->param;
334:   hx     = one / (PetscReal)(mx-1);
335:   hy     = one / (PetscReal)(my-1);
336:   sc     = hx*hy;
337:   hxdhy  = hx/hy;
338:   hydhx  = hy/hx;

340:   /*
341:      Get pointers to vector data
342:   */
343:   VecGetArray(X,&x);
344:   VecGetArray(F,&f);

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

367:   /*
368:      Restore vectors
369:   */
370:   VecRestoreArray(X,&x);
371:   VecRestoreArray(F,&f);
372:   return 0;
373: }
374: /* ------------------------------------------------------------------- */
377: /*
378:    FormJacobian - Evaluates Jacobian matrix.

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

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

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

408:   /*
409:      Get pointer to vector data
410:   */
411:   VecGetArray(X,&x);

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

432:   /*
433:      Restore vector
434:   */
435:   VecRestoreArray(X,&x);

437:   /* 
438:      Assemble matrix
439:   */
440:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
441:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

443:   if (jac != *J) {
444:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
445:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
446:   }

448:   /*
449:      Set flag to indicate that the Jacobian matrix retains an identical
450:      nonzero structure throughout all nonlinear iterations (although the
451:      values of the entries change). Thus, we can save some work in setting
452:      up the preconditioner (e.g., no need to redo symbolic factorization for
453:      ILU/ICC preconditioners).
454:       - If the nonzero structure of the matrix is different during
455:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
456:         must be used instead.  If you are unsure whether the matrix
457:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
458:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
459:         believes your assertion and does not check the structure
460:         of the matrix.  If you erroneously claim that the structure
461:         is the same when it actually is not, the new preconditioner
462:         will not function correctly.  Thus, use this optimization
463:         feature with caution!
464:   */
465:   *flag = SAME_NONZERO_PATTERN;
466:   return 0;
467: }