Actual source code: ex14.c
petsc-3.4.2 2013-07-02

      2: static char help[] = "Bratu nonlinear PDE in 3d.\n\
      3: We solve the  Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\
      4: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\
      5: The command line options include:\n\
      6:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
      7:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n";

      9: /*T
     10:    Concepts: SNES^parallel Bratu example
     11:    Concepts: DMDA^using distributed arrays;
     12:    Processors: n
     13: T*/

     15: /* ------------------------------------------------------------------------

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

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

     22:     with boundary conditions

     24:              u = 0  for  x = 0, x = 1, y = 0, y = 1, z = 0, z = 1

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

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

     33: /*
     34:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
     35:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
     36:    file automatically includes:
     37:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     38:      petscmat.h - matrices
     39:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     40:      petscviewer.h - viewers               petscpc.h  - preconditioners
     41:      petscksp.h   - linear solvers
     42: */
     43: #include <petscdmda.h>
     44: #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:   DM        da;                /* distributed array data structure */
     55: } AppCtx;

     57: /*
     58:    User-defined routines
     59: */
     60: extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*),FormInitialGuess(AppCtx*,Vec);
     61: extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);

     65: int main(int argc,char **argv)
     66: {
     67:   SNES           snes;                         /* nonlinear solver */
     68:   Vec            x,r;                          /* solution, residual vectors */
     69:   Mat            J;                            /* Jacobian matrix */
     70:   AppCtx         user;                         /* user-defined work context */
     71:   PetscInt       its;                          /* iterations for convergence */
     72:   MatFDColoring  matfdcoloring;
     73:   PetscBool      matrix_free = PETSC_FALSE,coloring = PETSC_FALSE;
     75:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm;

     77:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     78:      Initialize program
     79:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     81:   PetscInitialize(&argc,&argv,(char*)0,help);

     83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     84:      Initialize problem parameters
     85:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     86:   user.param = 6.0;
     87:   PetscOptionsGetReal(NULL,"-par",&user.param,NULL);
     88:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,1,"Lambda is out of range");

     90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     91:      Create nonlinear solver context
     92:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     93:   SNESCreate(PETSC_COMM_WORLD,&snes);

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:      Create distributed array (DMDA) to manage parallel grid and vectors
     97:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     98:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,-4,PETSC_DECIDE,PETSC_DECIDE,
     99:                       PETSC_DECIDE,1,1,NULL,NULL,NULL,&user.da);

    101:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    102:      Extract global vectors from DMDA; then duplicate for remaining
    103:      vectors that are the same types
    104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    105:   DMCreateGlobalVector(user.da,&x);
    106:   VecDuplicate(x,&r);

    108:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    109:      Set function evaluation routine and vector
    110:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    111:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

    113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    114:      Create matrix data structure; set Jacobian evaluation routine

    116:      Set Jacobian matrix data structure and default Jacobian evaluation
    117:      routine. User can override with:
    118:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
    119:                 (unless user explicitly sets preconditioner)
    120:      -snes_mf_operator : form preconditioning matrix as set by the user,
    121:                          but use matrix-free approx for Jacobian-vector
    122:                          products within Newton-Krylov method
    123:      -fdcoloring : using finite differences with coloring to compute the Jacobian

    125:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    126:   PetscOptionsGetBool(NULL,"-snes_mf",&matrix_free,NULL);
    127:   PetscOptionsGetBool(NULL,"-fdcoloring",&coloring,NULL);
    128:   if (!matrix_free) {
    129:     DMCreateMatrix(user.da,MATAIJ,&J);
    130:     if (coloring) {
    131:       ISColoring iscoloring;
    132:       DMCreateColoring(user.da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
    133:       MatFDColoringCreate(J,iscoloring,&matfdcoloring);
    134:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
    135:       MatFDColoringSetFromOptions(matfdcoloring);
    136:       SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
    137:       ISColoringDestroy(&iscoloring);
    138:     } else {
    139:       SNESSetJacobian(snes,J,J,FormJacobian,&user);
    140:     }
    141:   }

    143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    144:      Customize nonlinear solver; set runtime options
    145:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    146:   SNESSetFromOptions(snes);

    148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    149:      Evaluate initial guess
    150:      Note: The user should initialize the vector, x, with the initial guess
    151:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    152:      to employ an initial guess of zero, the user should explicitly set
    153:      this vector to zero by calling VecSet().
    154:   */
    155:   FormInitialGuess(&user,x);

    157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    158:      Solve nonlinear system
    159:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    160:   SNESSolve(snes,NULL,x);
    161:   SNESGetIterationNumber(snes,&its);

    163:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    164:      Explicitly check norm of the residual of the solution
    165:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    166:   FormFunction(snes,x,r,(void*)&user);
    167:   VecNorm(r,NORM_2,&fnorm);
    168:   PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %G\n",its,fnorm);

    170:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    171:      Free work space.  All PETSc objects should be destroyed when they
    172:      are no longer needed.
    173:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    175:   if (!matrix_free) {
    176:     MatDestroy(&J);
    177:   }
    178:   VecDestroy(&x);
    179:   VecDestroy(&r);
    180:   SNESDestroy(&snes);
    181:   DMDestroy(&user.da);
    182:   if (coloring) {MatFDColoringDestroy(&matfdcoloring);}
    183:   PetscFinalize();
    184:   return(0);
    185: }
    186: /* ------------------------------------------------------------------- */
    189: /*
    190:    FormInitialGuess - Forms initial approximation.

    192:    Input Parameters:
    193:    user - user-defined application context
    194:    X - vector

    196:    Output Parameter:
    197:    X - vector
    198:  */
    199: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
    200: {
    201:   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
    203:   PetscReal      lambda,temp1,hx,hy,hz,tempk,tempj;
    204:   PetscScalar    ***x;

    207:   DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    209:   lambda = user->param;
    210:   hx     = 1.0/(PetscReal)(Mx-1);
    211:   hy     = 1.0/(PetscReal)(My-1);
    212:   hz     = 1.0/(PetscReal)(Mz-1);
    213:   temp1  = lambda/(lambda + 1.0);

    215:   /*
    216:      Get a pointer to vector data.
    217:        - For default PETSc vectors, VecGetArray() returns a pointer to
    218:          the data array.  Otherwise, the routine is implementation dependent.
    219:        - You MUST call VecRestoreArray() when you no longer need access to
    220:          the array.
    221:   */
    222:   DMDAVecGetArray(user->da,X,&x);

    224:   /*
    225:      Get local grid boundaries (for 3-dimensional DMDA):
    226:        xs, ys, zs   - starting grid indices (no ghost points)
    227:        xm, ym, zm   - widths of local grid (no ghost points)

    229:   */
    230:   DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

    232:   /*
    233:      Compute initial guess over the locally owned part of the grid
    234:   */
    235:   for (k=zs; k<zs+zm; k++) {
    236:     tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz;
    237:     for (j=ys; j<ys+ym; j++) {
    238:       tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk);
    239:       for (i=xs; i<xs+xm; i++) {
    240:         if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
    241:           /* boundary conditions are all zero Dirichlet */
    242:           x[k][j][i] = 0.0;
    243:         } else {
    244:           x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj));
    245:         }
    246:       }
    247:     }
    248:   }

    250:   /*
    251:      Restore vector
    252:   */
    253:   DMDAVecRestoreArray(user->da,X,&x);
    254:   return(0);
    255: }
    256: /* ------------------------------------------------------------------- */
    259: /*
    260:    FormFunction - Evaluates nonlinear function, F(x).

    262:    Input Parameters:
    263: .  snes - the SNES context
    264: .  X - input vector
    265: .  ptr - optional user-defined context, as set by SNESSetFunction()

    267:    Output Parameter:
    268: .  F - function vector
    269:  */
    270: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
    271: {
    272:   AppCtx         *user = (AppCtx*)ptr;
    274:   PetscInt       i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm;
    275:   PetscReal      two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc;
    276:   PetscScalar    u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f;
    277:   Vec            localX;

    280:   DMGetLocalVector(user->da,&localX);
    281:   DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
    282:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    284:   lambda  = user->param;
    285:   hx      = 1.0/(PetscReal)(Mx-1);
    286:   hy      = 1.0/(PetscReal)(My-1);
    287:   hz      = 1.0/(PetscReal)(Mz-1);
    288:   sc      = hx*hy*hz*lambda;
    289:   hxhzdhy = hx*hz/hy;
    290:   hyhzdhx = hy*hz/hx;
    291:   hxhydhz = hx*hy/hz;

    293:   /*
    294:      Scatter ghost points to local vector,using the 2-step process
    295:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    296:      By placing code between these two statements, computations can be
    297:      done while messages are in transition.
    298:   */
    299:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
    300:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

    302:   /*
    303:      Get pointers to vector data
    304:   */
    305:   DMDAVecGetArray(user->da,localX,&x);
    306:   DMDAVecGetArray(user->da,F,&f);

    308:   /*
    309:      Get local grid boundaries
    310:   */
    311:   DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

    313:   /*
    314:      Compute function over the locally owned part of the grid
    315:   */
    316:   for (k=zs; k<zs+zm; k++) {
    317:     for (j=ys; j<ys+ym; j++) {
    318:       for (i=xs; i<xs+xm; i++) {
    319:         if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) {
    320:           f[k][j][i] = x[k][j][i];
    321:         } else {
    322:           u          = x[k][j][i];
    323:           u_east     = x[k][j][i+1];
    324:           u_west     = x[k][j][i-1];
    325:           u_north    = x[k][j+1][i];
    326:           u_south    = x[k][j-1][i];
    327:           u_up       = x[k+1][j][i];
    328:           u_down     = x[k-1][j][i];
    329:           u_xx       = (-u_east + two*u - u_west)*hyhzdhx;
    330:           u_yy       = (-u_north + two*u - u_south)*hxhzdhy;
    331:           u_zz       = (-u_up + two*u - u_down)*hxhydhz;
    332:           f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u);
    333:         }
    334:       }
    335:     }
    336:   }

    338:   /*
    339:      Restore vectors
    340:   */
    341:   DMDAVecRestoreArray(user->da,localX,&x);
    342:   DMDAVecRestoreArray(user->da,F,&f);
    343:   DMRestoreLocalVector(user->da,&localX);
    344:   PetscLogFlops(11.0*ym*xm);
    345:   return(0);
    346: }
    347: /* ------------------------------------------------------------------- */
    350: /*
    351:    FormJacobian - Evaluates Jacobian matrix.

    353:    Input Parameters:
    354: .  snes - the SNES context
    355: .  x - input vector
    356: .  ptr - optional user-defined context, as set by SNESSetJacobian()

    358:    Output Parameters:
    359: .  A - Jacobian matrix
    360: .  B - optionally different preconditioning matrix
    361: .  flag - flag indicating matrix structure

    363: */
    364: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
    365: {
    366:   AppCtx         *user = (AppCtx*)ptr;  /* user-defined application context */
    367:   Mat            jac   = *B;              /* Jacobian matrix */
    368:   Vec            localX;
    370:   PetscInt       i,j,k,Mx,My,Mz;
    371:   MatStencil     col[7],row;
    372:   PetscInt       xs,ys,zs,xm,ym,zm;
    373:   PetscScalar    lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x;

    376:   DMGetLocalVector(user->da,&localX);
    377:   DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,
    378:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    380:   lambda  = user->param;
    381:   hx      = 1.0/(PetscReal)(Mx-1);
    382:   hy      = 1.0/(PetscReal)(My-1);
    383:   hz      = 1.0/(PetscReal)(Mz-1);
    384:   sc      = hx*hy*hz*lambda;
    385:   hxhzdhy = hx*hz/hy;
    386:   hyhzdhx = hy*hz/hx;
    387:   hxhydhz = hx*hy/hz;

    389:   /*
    390:      Scatter ghost points to local vector, using the 2-step process
    391:         DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
    392:      By placing code between these two statements, computations can be
    393:      done while messages are in transition.
    394:   */
    395:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
    396:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);

    398:   /*
    399:      Get pointer to vector data
    400:   */
    401:   DMDAVecGetArray(user->da,localX,&x);

    403:   /*
    404:      Get local grid boundaries
    405:   */
    406:   DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);

    408:   /*
    409:      Compute entries for the locally owned part of the Jacobian.
    410:       - Currently, all PETSc parallel matrix formats are partitioned by
    411:         contiguous chunks of rows across the processors.
    412:       - Each processor needs to insert only elements that it owns
    413:         locally (but any non-local elements will be sent to the
    414:         appropriate processor during matrix assembly).
    415:       - Here, we set all entries for a particular row at once.
    416:       - We can set matrix entries either using either
    417:         MatSetValuesLocal() or MatSetValues(), as discussed above.
    418:   */
    419:   for (k=zs; k<zs+zm; k++) {
    420:     for (j=ys; j<ys+ym; j++) {
    421:       for (i=xs; i<xs+xm; i++) {
    422:         row.k = k; row.j = j; row.i = i;
    423:         /* boundary points */
    424:         if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) {
    425:           v[0] = 1.0;
    426:           MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
    427:         } else {
    428:           /* interior grid points */
    429:           v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j;  col[0].i = i;
    430:           v[1] = -hxhzdhy; col[1].k=k;  col[1].j=j-1;col[1].i = i;
    431:           v[2] = -hyhzdhx; col[2].k=k;  col[2].j=j;  col[2].i = i-1;
    432:           v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i;
    433:           v[4] = -hyhzdhx; col[4].k=k;  col[4].j=j;  col[4].i = i+1;
    434:           v[5] = -hxhzdhy; col[5].k=k;  col[5].j=j+1;col[5].i = i;
    435:           v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j;  col[6].i = i;
    436:           MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
    437:         }
    438:       }
    439:     }
    440:   }
    441:   DMDAVecRestoreArray(user->da,localX,&x);
    442:   DMRestoreLocalVector(user->da,&localX);

    444:   /*
    445:      Assemble matrix, using the 2-step process:
    446:        MatAssemblyBegin(), MatAssemblyEnd().
    447:   */
    448:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    449:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

    451:   /*
    452:      Normally since the matrix has already been assembled above; this
    453:      would do nothing. But in the matrix free mode -snes_mf_operator
    454:      this tells the "matrix-free" matrix that a new linear system solve
    455:      is about to be done.
    456:   */

    458:   MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
    459:   MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);

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

    481:   /*
    482:      Tell the matrix we will never add a new nonzero location to the
    483:      matrix. If we do, it will generate an error.
    484:   */
    485:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
    486:   return(0);
    487: }