Actual source code: ex1.c

petsc-3.11.4 2019-09-28
Report Typos and Errors

  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*/



 16: /* ------------------------------------------------------------------------

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

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

 23:     with boundary conditions

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

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

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

 33:   ------------------------------------------------------------------------- */

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

 45:  #include <petscsnes.h>

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

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

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

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

 86:   /*
 87:      Initialize problem parameters
 88:   */
 89:   user.mx = 4; user.my = 4; user.param = 6.0;
 90:   PetscOptionsGetInt(NULL,NULL,"-mx",&user.mx,NULL);
 91:   PetscOptionsGetInt(NULL,NULL,"-my",&user.my,NULL);
 92:   PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);
 93:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"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);



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) {
140:     MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5,NULL,&J);
141:   }

143:   /*
144:      This option will cause the Jacobian to be computed via finite differences
145:     efficiently using a coloring of the columns of the matrix.
146:   */
147:   PetscOptionsGetBool(NULL,NULL,"-snes_fd_coloring",&fd_coloring,NULL);
148:   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");

150:   if (fd_coloring) {
151:     ISColoring   iscoloring;
152:     MatColoring  mc;

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

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

204:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205:      Customize nonlinear solver; set runtime options
206:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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


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

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

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

278:    Input Parameters:
279:    user - user-defined application context
280:    X - vector

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

292:   mx     = user->mx;
293:   my     = user->my;
294:   lambda = user->param;

296:   hx = 1.0 / (PetscReal)(mx-1);
297:   hy = 1.0 / (PetscReal)(my-1);

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

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

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

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

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

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

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

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

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

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

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

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

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

446:   /*
447:      Restore vector
448:   */
449:   VecRestoreArrayRead(X,&x);

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

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

462:   return 0;
463: }

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

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

478: PetscErrorCode ConvergenceDestroy(void* ctx)
479: {

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



490: /*TEST

492:    build:
493:       requires: !single

495:    test:
496:       args: -ksp_gmres_cgs_refinement_type refine_always

498:    test:
499:       suffix: 2
500:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always

502:    test:
503:       suffix: 2a
504:       filter: grep -i KSPConvergedDefault > /dev/null && echo "Found KSPConvergedDefault"
505:       args: -snes_monitor_short -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always -info
506:       requires: define(PETSC_USE_LOG)

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

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

518: TEST*/