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

      1: static char help[] = "Solves the Lane-Emden equation in a 2D rectangular\n\
      2: domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\n";

      4: /*T
      5:    Concepts: SNES^parallel Lane-Emden example
      6:    Concepts: DMDA^using distributed arrays;
      7:    Processors: n
      8: T*/

     10: /* ------------------------------------------------------------------------

     12:     The Lane-Emden equation is given by the partial differential equation

     14:             -alpha*Laplacian u - lambda*u^3 = 0,  0 < x,y < 1,

     16:     with boundary conditions

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

     20:     A bilinear finite element approximation is used to discretize the boundary
     21:     value problem to obtain a nonlinear system of equations.

     23:   ------------------------------------------------------------------------- */

     25: /*
     26:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
     27:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
     28:    file automatically includes:
     29:      petscsys.h       - base PETSc routines   petscvec.h - vectors
     30:      petscmat.h - matrices
     31:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     32:      petscviewer.h - viewers               petscpc.h  - preconditioners
     33:      petscksp.h   - linear solvers
     34: */
     35: #include <petscsys.h>
     36: #include <petscbag.h>
     37: #include <petscdmda.h>
     38: #include <petscsnes.h>

     40: /*
     41:    User-defined application context - contains data needed by the
     42:    application-provided call-back routines, FormJacobianLocal() and
     43:    FormFunctionLocal().
     44: */
     45: typedef struct {
     46:   PetscReal alpha;           /* parameter controlling linearity */
     47:   PetscReal lambda;          /* parameter controlling nonlinearity */
     48: } AppCtx;

     50: static PetscScalar Kref[16] = { 0.666667, -0.166667, -0.333333, -0.166667,
     51:                                -0.166667,  0.666667, -0.166667, -0.333333,
     52:                                -0.333333, -0.166667,  0.666667, -0.166667,
     53:                                -0.166667, -0.333333, -0.166667,  0.666667};

     55: /* These are */
     56: static PetscScalar quadPoints[8] = {0.211325, 0.211325,
     57:                                     0.788675, 0.211325,
     58:                                     0.788675, 0.788675,
     59:                                     0.211325, 0.788675};
     60: static PetscScalar quadWeights[4] = {0.25, 0.25, 0.25, 0.25};

     62: /*
     63:    User-defined routines
     64: */
     65: extern PetscErrorCode FormInitialGuess(SNES,Vec,void*);
     66: extern PetscErrorCode FormFunctionLocal(DMDALocalInfo*,PetscScalar**,PetscScalar**,AppCtx*);
     67: extern PetscErrorCode FormJacobianLocal(DMDALocalInfo*,PetscScalar**,Mat,Mat,MatStructure*,AppCtx*);
     68: extern PetscErrorCode PrintVector(DM, Vec);

     72: int main(int argc,char **argv)
     73: {
     74:   DM                  da;
     75:   SNES                snes;                    /* nonlinear solver */
     76:   AppCtx              *user;                   /* user-defined work context */
     77:   PetscBag            bag;
     78:   PetscInt            its;                     /* iterations for convergence */
     79:   SNESConvergedReason reason;
     80:   PetscErrorCode      ierr;
     81:   PetscReal           lambda_max = 6.81, lambda_min = 0.0;
     82:   Vec                 x;

     84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     85:      Initialize program
     86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     87:   PetscInitialize(&argc,&argv,(char*)0,help);

     89:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     90:      Initialize problem parameters
     91:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     92:   PetscBagCreate(PETSC_COMM_WORLD, sizeof(AppCtx), &bag);
     93:   PetscBagGetData(bag, (void**) &user);
     94:   PetscBagSetName(bag, "params", "Parameters for SNES example 4");
     95:   PetscBagRegisterReal(bag, &user->alpha, 1.0, "alpha", "Linear coefficient");
     96:   PetscBagRegisterReal(bag, &user->lambda, 6.0, "lambda", "Nonlinear coefficient");
     97:   PetscBagSetFromOptions(bag);
     98:   PetscOptionsGetReal(NULL,"-alpha",&user->alpha,NULL);
     99:   PetscOptionsGetReal(NULL,"-lambda",&user->lambda,NULL);
    100:   if (user->lambda > lambda_max || user->lambda < lambda_min) SETERRQ3(PETSC_COMM_SELF,1,"Lambda %G is out of range [%G, %G]", user->lambda, lambda_min, lambda_max);

    102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    103:      Create SNES to manage hierarchical solvers
    104:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    105:   SNESCreate(PETSC_COMM_WORLD,&snes);

    107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    108:      Create distributed array (DMDA) to manage parallel grid and vectors
    109:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    110:   DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,-3,-3,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
    111:   DMDASetFieldName(da, 0, "ooblek");
    112:   DMSetApplicationContext(da,user);
    113:   SNESSetDM(snes, (DM) da);

    115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    116:      Set the discretization functions
    117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    118:   DMDASNESSetFunctionLocal(da,INSERT_VALUES,(PetscErrorCode (*)(DMDALocalInfo*,void*,void*,void*))FormFunctionLocal,user);
    119:   DMDASNESSetJacobianLocal(da,(PetscErrorCode (*)(DMDALocalInfo*,void*,Mat,Mat,MatStructure*,void*))FormJacobianLocal,user);
    120:   SNESSetFromOptions(snes);

    122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    123:      Evaluate initial guess
    124:      Note: The user should initialize the vector, x, with the initial guess
    125:      for the nonlinear solver prior to calling SNESSolve().  In particular,
    126:      to employ an initial guess of zero, the user should explicitly set
    127:      this vector to zero by calling VecSet().
    128:   */
    129:   SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL);

    131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    132:      Solve nonlinear system
    133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    134:   SNESSolve(snes,NULL,NULL);
    135:   SNESGetIterationNumber(snes,&its);
    136:   SNESGetConvergedReason(snes, &reason);
    137:   DMDestroy(&da);
    138:   SNESGetDM(snes,&da);
    139:   SNESGetSolution(snes,&x);
    140:   PrintVector(da, x);

    142:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    143:      Free work space.  All PETSc objects should be destroyed when they
    144:      are no longer needed.
    145:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    146:   SNESDestroy(&snes);
    147:   PetscBagDestroy(&bag);
    148:   PetscFinalize();
    149:   return(0);
    150: }

    154: PetscErrorCode PrintVector(DM da, Vec U)
    155: {
    156:   PetscScalar    **u;
    157:   PetscInt       i,j,xs,ys,xm,ym;

    161:   DMDAVecGetArray(da,U,&u);
    162:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
    163:   for (j = ys+ym-1; j >= ys; j--) {
    164:     for (i = xs; i < xs+xm; i++) {
    165:       PetscPrintf(PETSC_COMM_SELF,"u[%d][%d] = %G ", j, i, PetscRealPart(u[j][i]));
    166:     }
    167:     PetscPrintf(PETSC_COMM_SELF,"\n");
    168:   }
    169:   DMDAVecRestoreArray(da,U,&u);
    170:   return(0);
    171: }

    175: PetscErrorCode ExactSolution(PetscReal x, PetscReal y, PetscScalar *u)
    176: {
    178:   *u = x*x;
    179:   return(0);
    180: }

    184: /*
    185:    FormInitialGuess - Forms initial approximation.

    187:    Input Parameters:
    188:    X - vector

    190:    Output Parameter:
    191:    X - vector
    192: */
    193: PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx)
    194: {
    195:   AppCtx         *user;
    196:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
    198:   PetscReal      lambda,temp1,temp,hx,hy;
    199:   PetscScalar    **x;
    200:   DM             da;

    203:   SNESGetDM(snes,&da);
    204:   DMGetApplicationContext(da,&user);
    205:   DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
    206:                      PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

    208:   lambda = user->lambda;
    209:   hx     = 1.0/(PetscReal)(Mx-1);
    210:   hy     = 1.0/(PetscReal)(My-1);
    211:   if (lambda == 0.0) temp1 = 1.0;
    212:   else temp1 = lambda/(lambda + 1.0);

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

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

    228:   */
    229:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

    231:   /*
    232:      Compute initial guess over the locally owned part of the grid
    233:   */
    234:   for (j=ys; j<ys+ym; j++) {
    235:     temp = (PetscReal)(PetscMin(j,My-j-1))*hy;
    236:     for (i=xs; i<xs+xm; i++) {

    238:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) x[j][i] = 0.0; /* boundary conditions are all zero Dirichlet */
    239:       else x[j][i] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,temp));
    240:     }
    241:   }

    243:   /*
    244:      Restore vector
    245:   */
    246:   DMDAVecRestoreArray(da,X,&x);
    247:   return(0);
    248: }

    252: PetscErrorCode constantResidual(PetscReal lambda, int i, int j, PetscReal hx, PetscReal hy, PetscScalar r[])
    253: {
    254:   PetscScalar rLocal[4] = {0.0, 0.0, 0.0};
    255:   PetscScalar phi[4]    = {0.0, 0.0, 0.0, 0.0};
    256:   /* PetscReal   xI = i*hx, yI = j*hy, x, y; */
    257:   PetscReal   hxhy = hx*hy;
    258:   PetscScalar res;
    259:   PetscInt    q, k;

    262:   for (q = 0; q < 4; q++) {
    263:     phi[0] = (1.0 - quadPoints[q*2])*(1.0 - quadPoints[q*2+1]);
    264:     phi[1] =  quadPoints[q*2]       *(1.0 - quadPoints[q*2+1]);
    265:     phi[2] =  quadPoints[q*2]       * quadPoints[q*2+1];
    266:     phi[3] = (1.0 - quadPoints[q*2])* quadPoints[q*2+1];
    267:     /*
    268:      x      = xI + quadPoints[q*2]*hx;
    269:      y      = yI + quadPoints[q*2+1]*hy;
    270:      */
    271:     res = quadWeights[q]*(2.0);
    272:     for (k = 0; k < 4; k++) rLocal[k] += phi[k]*res;
    273:   }
    274:   for (k = 0; k < 4; k++) r[k] += lambda*hxhy*rLocal[k];
    275:   return(0);
    276: }

    280: PetscErrorCode nonlinearResidual(PetscReal lambda, PetscScalar u[], PetscScalar r[])
    281: {
    283:   r[0] += lambda*(48.0*u[0]*u[0]*u[0] + 12.0*u[1]*u[1]*u[1] + 9.0*u[0]*u[0]*(4.0*u[1] + u[2] + 4.0*u[3]) + u[1]*u[1]*(9.0*u[2] + 6.0*u[3]) + u[1]*(6.0*u[2]*u[2] + 8.0*u[2]*u[3] + 6.0*u[3]*u[3])
    284:                   + 3.0*(u[2]*u[2]*u[2] + 2.0*u[2]*u[2]*u[3] + 3.0*u[2]*u[3]*u[3] + 4.0*u[3]*u[3]*u[3])
    285:                   + 2.0*u[0]*(12.0*u[1]*u[1] + u[1]*(6.0*u[2] + 9.0*u[3]) + 2.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3])))/1200.0;
    286:   r[1] += lambda*(12.0*u[0]*u[0]*u[0] + 48.0*u[1]*u[1]*u[1] + 9.0*u[1]*u[1]*(4.0*u[2] + u[3]) + 3.0*u[0]*u[0]*(8.0*u[1] + 2.0*u[2] + 3.0*u[3])
    287:                   + 4.0*u[1]*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3]) + 3.0*(4.0*u[2]*u[2]*u[2] + 3.0*u[2]*u[2]*u[3] + 2.0*u[2]*u[3]*u[3] + u[3]*u[3]*u[3])
    288:                   + 2.0*u[0]*((18.0*u[1]*u[1] + 3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3]) + u[1]*(9.0*u[2] + 6.0*u[3])))/1200.0;
    289:   r[2] += lambda*(3.0*u[0]*u[0]*u[0] + u[0]*u[0]*(6.0*u[1] + 4.0*u[2] + 6.0*u[3]) + u[0]*(9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]))
    290:                   + 6.0*(2.0*u[1]*u[1]*u[1] + u[1]*u[1]*(4.0*u[2] + u[3]) + u[1]*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3]) + 2.0*(4.0*u[2]*u[2]*u[2] + 3.0*u[2]*u[2]*u[3] + 2.0*u[2]*u[3]*u[3] + u[3]*u[3]*u[3])))/1200.0;
    291:   r[3] += lambda*(12.0*u[0]*u[0]*u[0] + 3.0*u[1]*u[1]*u[1] + u[1]*u[1]*(6.0*u[2] + 4.0*u[3]) + 3.0*u[0]*u[0]*(3.0*u[1] + 2.0*u[2] + 8.0*u[3])
    292:                   + 3.0*u[1]*(3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3]) + 12.0*(u[2]*u[2]*u[2] + 2.0*u[2]*u[2]*u[3] + 3.0*u[2]*u[3]*u[3] + 4.0*u[3]*u[3]*u[3])
    293:                   + 2.0*u[0]*(3.0*u[1]*u[1] + u[1]*(4.0*u[2] + 6.0*u[3]) + 3.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3])))/1200.0;
    294:   return(0);
    295: }

    299: PetscErrorCode nonlinearResidualBratu(PetscReal lambda, PetscScalar u[], PetscScalar r[])
    300: {
    301:   PetscScalar rLocal[4] = {0.0, 0.0, 0.0, 0.0};
    302:   PetscScalar phi[4]    = {0.0, 0.0, 0.0, 0.0};
    303:   PetscScalar res;
    304:   PetscInt    q;

    307:   for (q = 0; q < 4; q++) {
    308:     phi[0]     = (1.0 - quadPoints[q*2])*(1.0 - quadPoints[q*2+1]);
    309:     phi[1]     =  quadPoints[q*2]       *(1.0 - quadPoints[q*2+1]);
    310:     phi[2]     =  quadPoints[q*2]       * quadPoints[q*2+1];
    311:     phi[3]     = (1.0 - quadPoints[q*2])* quadPoints[q*2+1];
    312:     res        = quadWeights[q]*PetscExpScalar(u[0]*phi[0]+ u[1]*phi[1] + u[2]*phi[2]+ u[3]*phi[3]);
    313:     rLocal[0] += phi[0]*res;
    314:     rLocal[1] += phi[1]*res;
    315:     rLocal[2] += phi[2]*res;
    316:     rLocal[3] += phi[3]*res;
    317:   }
    318:   r[0] += lambda*rLocal[0];
    319:   r[1] += lambda*rLocal[1];
    320:   r[2] += lambda*rLocal[2];
    321:   r[3] += lambda*rLocal[3];
    322:   return(0);
    323: }

    327: PetscErrorCode nonlinearJacobian(PetscScalar lambda, PetscScalar u[], PetscScalar J[])
    328: {
    330:   J[0]  = lambda*(72.0*u[0]*u[0] + 12.0*u[1]*u[1] + 9.0*u[0]*(4.0*u[1] + u[2] + 4.0*u[3]) + u[1]*(6.0*u[2] + 9.0*u[3]) + 2.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]))/600.0;
    331:   J[1]  = lambda*(18.0*u[0]*u[0] + 18.0*u[1]*u[1] + 3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3] + 3.0*u[0]*(8.0*u[1] + 2.0*u[2] + 3.0*u[3]) + u[1]*(9.0*u[2] + 6.0*u[3]))/600.0;
    332:   J[2]  = lambda*( 9.0*u[0]*u[0] +  9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
    333:   J[3]  = lambda*(18.0*u[0]*u[0] +  3.0*u[1]*u[1] + u[1]*(4.0*u[2] + 6.0*u[3]) + 3.0*u[0]*(3.0*u[1] + 2.0*u[2] + 8.0*u[3]) + 3.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]))/600.0;

    335:   J[4]  = lambda*(18.0*u[0]*u[0] + 18.0*u[1]*u[1] + 3.0*u[2]*u[2] + 4.0*u[2]*u[3] + 3.0*u[3]*u[3] + 3.0*u[0]*(8.0*u[1] + 2.0*u[2] + 3.0*u[3]) + u[1]*(9.0*u[2] + 6.0*u[3]))/600.0;
    336:   J[5]  = lambda*(12.0*u[0]*u[0] + 72.0*u[1]*u[1] + 9.0*u[1]*(4.0*u[2] + u[3]) + u[0]*(36.0*u[1] + 9.0*u[2] + 6.0*u[3]) + 2.0*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3]))/600.0;
    337:   J[6]  = lambda*( 3.0*u[0]*u[0] + u[0]*(9.0*u[1] + 6.0*u[2] + 4.0*u[3]) + 3.0*(6.0*u[1]*u[1] + 6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3] + 2.0*u[1]*(4.0*u[2] + u[3])))/600.0;
    338:   J[7]  = lambda*( 9.0*u[0]*u[0] +  9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;

    340:   J[8]  = lambda*( 9.0*u[0]*u[0] +  9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
    341:   J[9]  = lambda*( 3.0*u[0]*u[0] + u[0]*(9.0*u[1] + 6.0*u[2] + 4.0*u[3]) + 3.0*(6.0*u[1]*u[1] + 6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3] + 2.0*u[1]*(4.0*u[2] + u[3])))/600.0;
    342:   J[10] = lambda*( 2.0*u[0]*u[0] + u[0]*(6.0*u[1] + 9.0*u[2] + 6.0*u[3]) + 3.0*(4.0*u[1]*u[1] + 3.0*u[1]*(4.0*u[2] + u[3]) + 4.0*(6.0*u[2]*u[2] + 3.0*u[2]*u[3] + u[3]*u[3])))/600.0;
    343:   J[11] = lambda*( 3.0*u[0]*u[0] + u[0]*(4.0*u[1] + 6.0*u[2] + 9.0*u[3]) + 3.0*(u[1]*u[1] + 6.0*u[2]*u[2] + 8.0*u[2]*u[3] + 6.0*u[3]*u[3] + u[1]*(3.0*u[2] + 2.0*u[3])))/600.0;

    345:   J[12] = lambda*(18.0*u[0]*u[0] +  3.0*u[1]*u[1] + u[1]*(4.0*u[2] + 6.0*u[3]) + 3.0*u[0]*(3.0*u[1] + 2.0*u[2] + 8.0*u[3]) + 3.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]))/600.0;
    346:   J[13] = lambda*( 9.0*u[0]*u[0] +  9.0*u[1]*u[1] + 9.0*u[2]*u[2] + 12.0*u[2]*u[3] + 9.0*u[3]*u[3] + 4.0*u[1]*(3.0*u[2] + 2.0*u[3]) + 4.0*u[0]*(3.0*u[1] + 2.0*u[2] + 3.0*u[3]))/1200.0;
    347:   J[14] = lambda*( 3.0*u[0]*u[0] + u[0]*(4.0*u[1] + 6.0*u[2] + 9.0*u[3]) + 3.0*(u[1]*u[1] + 6.0*u[2]*u[2] + 8.0*u[2]*u[3] + 6.0*u[3]*u[3] + u[1]*(3.0*u[2] + 2.0*u[3])))/600.0;
    348:   J[15] = lambda*(12.0*u[0]*u[0] +  2.0*u[1]*u[1] + u[1]*(6.0*u[2] + 9.0*u[3]) + 12.0*(u[2]*u[2] + 3.0*u[2]*u[3] + 6.0*u[3]*u[3]) + u[0]*(6.0*u[1] + 9.0*(u[2] + 4.0*u[3])))/600.0;
    349:   return(0);
    350: }

    354: /*
    355:    FormFunctionLocal - Evaluates nonlinear function, F(x).

    357:  */
    358: PetscErrorCode FormFunctionLocal(DMDALocalInfo *info,PetscScalar **x,PetscScalar **f,AppCtx *user)
    359: {
    360:   PetscScalar    uLocal[4];
    361:   PetscScalar    rLocal[4];
    362:   PetscScalar    uExact;
    363:   PetscReal      alpha,lambda,hx,hy,hxhy,sc;
    364:   PetscInt       i,j,k,l;

    368:   alpha  = user->alpha;
    369:   lambda = user->lambda;
    370:   hx     = 1.0/(PetscReal)(info->mx-1);
    371:   hy     = 1.0/(PetscReal)(info->my-1);
    372:   sc     = hx*hy*lambda;
    373:   hxhy   = hx*hy;

    375:   /* Zero the vector */
    376:   PetscMemzero((void*) &(f[info->xs][info->ys]), info->xm*info->ym*sizeof(PetscScalar));
    377:   /* Compute function over the locally owned part of the grid. For each
    378:      vertex (i,j), we consider the element below:

    380:        3         2
    381:      i,j+1 --- i+1,j+1
    382:        |         |
    383:        |         |
    384:       i,j  --- i+1,j
    385:        0         1

    387:      and therefore we do not loop over the last vertex in each dimension.
    388:   */
    389:   for (j = info->ys; j < info->ys+info->ym-1; j++) {
    390:     for (i = info->xs; i < info->xs+info->xm-1; i++) {
    391:       uLocal[0] = x[j][i];
    392:       uLocal[1] = x[j][i+1];
    393:       uLocal[2] = x[j+1][i+1];
    394:       uLocal[3] = x[j+1][i];
    395:       for (k = 0; k < 4; k++) {
    396:         rLocal[k] = 0.0;
    397:         for (l = 0; l < 4; l++) rLocal[k] += Kref[k*4 + l]*uLocal[l];
    398:         rLocal[k] *= hxhy*alpha;
    399:       }
    400:       constantResidual(1.0, i, j, hx, hy, rLocal);
    401:       nonlinearResidual(-1.0*sc, uLocal, rLocal);

    403:       f[j][i]     += rLocal[0];
    404:       f[j][i+1]   += rLocal[1];
    405:       f[j+1][i+1] += rLocal[2];
    406:       f[j+1][i]   += rLocal[3];

    408:       if (i == 0 || j == 0) {
    409:         ExactSolution(i*hx, j*hy, &uExact);

    411:         f[j][i] = x[j][i] - uExact;
    412:       }
    413:       if ((i == info->mx-2) || (j == 0)) {
    414:         ExactSolution((i+1)*hx, j*hy, &uExact);

    416:         f[j][i+1] = x[j][i+1] - uExact;
    417:       }
    418:       if ((i == info->mx-2) || (j == info->my-2)) {
    419:         ExactSolution((i+1)*hx, (j+1)*hy, &uExact);

    421:         f[j+1][i+1] = x[j+1][i+1] - uExact;
    422:       }
    423:       if ((i == 0) || (j == info->my-2)) {
    424:         ExactSolution(i*hx, (j+1)*hy, &uExact);

    426:         f[j+1][i] = x[j+1][i] - uExact;
    427:       }
    428:     }
    429:   }

    431:   PetscLogFlops(68.0*(info->ym-1)*(info->xm-1));
    432:   return(0);
    433: }

    437: /*
    438:    FormJacobianLocal - Evaluates Jacobian matrix.
    439: */
    440: PetscErrorCode FormJacobianLocal(DMDALocalInfo *info,PetscScalar **x,Mat A,Mat jac,MatStructure *str,AppCtx *user)
    441: {
    442:   PetscScalar    JLocal[16], ELocal[16], uLocal[4];
    443:   MatStencil     rows[4], cols[4], ident;
    444:   PetscInt       localRows[4];
    445:   PetscScalar    alpha,lambda,hx,hy,hxhy,sc;
    446:   PetscInt       i,j,k,l,numRows;

    450:   alpha  = user->alpha;
    451:   lambda = user->lambda;
    452:   hx     = 1.0/(PetscReal)(info->mx-1);
    453:   hy     = 1.0/(PetscReal)(info->my-1);
    454:   sc     = hx*hy*lambda;
    455:   hxhy   = hx*hy;

    457:   MatZeroEntries(jac);
    458:   /*
    459:      Compute entries for the locally owned part of the Jacobian.
    460:       - Currently, all PETSc parallel matrix formats are partitioned by
    461:         contiguous chunks of rows across the processors.
    462:       - Each processor needs to insert only elements that it owns
    463:         locally (but any non-local elements will be sent to the
    464:         appropriate processor during matrix assembly).
    465:       - Here, we set all entries for a particular row at once.
    466:       - We can set matrix entries either using either
    467:         MatSetValuesLocal() or MatSetValues(), as discussed above.
    468:   */
    469:   for (j=info->ys; j<info->ys+info->ym-1; j++) {
    470:     for (i=info->xs; i<info->xs+info->xm-1; i++) {
    471:       numRows   = 0;
    472:       uLocal[0] = x[j][i];
    473:       uLocal[1] = x[j][i+1];
    474:       uLocal[2] = x[j+1][i+1];
    475:       uLocal[3] = x[j+1][i];
    476:       /* i,j */
    477:       if (i == 0 || j == 0) {
    478:         ident.i   = i; ident.j = j;
    479:         JLocal[0] = 1.0;

    481:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    482:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    483:         MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
    484:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    485:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    486:       } else {
    487:         localRows[numRows] = 0;
    488:         rows[numRows].i    = i; rows[numRows].j = j;
    489:         numRows++;
    490:       }
    491:       cols[0].i = i; cols[0].j = j;
    492:       /* i+1,j */
    493:       if ((i == info->mx-2) || (j == 0)) {
    494:         ident.i   = i+1; ident.j = j;
    495:         JLocal[0] = 1.0;

    497:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    498:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    499:         MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
    500:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    501:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    502:       } else {
    503:         localRows[numRows] = 1;
    504:         rows[numRows].i    = i+1; rows[numRows].j = j;
    505:         numRows++;
    506:       }
    507:       cols[1].i = i+1; cols[1].j = j;
    508:       /* i+1,j+1 */
    509:       if ((i == info->mx-2) || (j == info->my-2)) {
    510:         ident.i   = i+1; ident.j = j+1;
    511:         JLocal[0] = 1.0;

    513:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    514:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    515:         MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
    516:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    517:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    518:       } else {
    519:         localRows[numRows] = 2;
    520:         rows[numRows].i    = i+1; rows[numRows].j = j+1;
    521:         numRows++;
    522:       }
    523:       cols[2].i = i+1; cols[2].j = j+1;
    524:       /* i,j+1 */
    525:       if ((i == 0) || (j == info->my-2)) {
    526:         ident.i   = i; ident.j = j+1;
    527:         JLocal[0] = 1.0;

    529:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    530:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    531:         MatSetValuesStencil(jac,1,&ident,1,&ident,JLocal,INSERT_VALUES);
    532:         MatAssemblyBegin(jac,MAT_FLUSH_ASSEMBLY);
    533:         MatAssemblyEnd(jac,MAT_FLUSH_ASSEMBLY);
    534:       } else {
    535:         localRows[numRows] = 3;
    536:         rows[numRows].i    = i; rows[numRows].j = j+1;
    537:         numRows++;
    538:       }
    539:       cols[3].i = i; cols[3].j = j+1;
    540:       nonlinearJacobian(-1.0*sc, uLocal, ELocal);
    541:       for (k = 0; k < numRows; k++) {
    542:         for (l = 0; l < 4; l++) {
    543:           JLocal[k*4 + l] = (hxhy*alpha*Kref[localRows[k]*4 + l] + ELocal[localRows[k]*4 + l]);
    544:         }
    545:       }
    546:       MatSetValuesStencil(jac,numRows,rows,4,cols,JLocal,ADD_VALUES);
    547:     }
    548:   }

    550:   /*
    551:      Assemble matrix, using the 2-step process:
    552:        MatAssemblyBegin(), MatAssemblyEnd().
    553:   */
    554:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    555:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
    556:   if (A != jac) {
    557:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
    558:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
    559:   }
    560:   *str = SAME_NONZERO_PATTERN;

    562:   /*
    563:      Tell the matrix we will never add a new nonzero location to the
    564:      matrix. If we do, it will generate an error.
    565:   */
    566:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
    567:   return(0);
    568: }