Actual source code: ex13.c

  2: static char help[] = "This program is a replica of ex6.c except that it does 2 solves to avoid paging.\n\
  3: This program demonstrates use of the SNES package to solve systems of\n\
  4: nonlinear equations in parallel, using 2-dimensional distributed arrays.\n\
  5: The 2-dim Bratu (SFI - solid fuel ignition) test problem is used, where\n\
  6: analytic formation of the Jacobian is the default.  The command line\n\
  7: options are:\n\
  8:   -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  9:      problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
 10:   -mx <xg>, where <xg> = number of grid points in the x-direction\n\
 11:   -my <yg>, where <yg> = number of grid points in the y-direction\n\
 12:   -Nx <npx>, where <npx> = number of processors in the x-direction\n\
 13:   -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";

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

 30:  #include petscsnes.h
 31:  #include petscda.h

 33: /* User-defined application context */
 34: typedef struct {
 35:    PetscReal   param;         /* test problem parameter */
 36:    PetscInt    mx,my;         /* discretization in x, y directions */
 37:    Vec         localX,localF; /* ghosted local vector */
 38:    DA          da;            /* distributed array data structure */
 39: } AppCtx;


 46: int main(int argc,char **argv)
 47: {
 48:   SNES           snes;                      /* nonlinear solver */
 49:   const SNESType type = SNESLS;             /* nonlinear solution method */
 50:   Vec            x,r;                       /* solution, residual vectors */
 51:   Mat            J;                         /* Jacobian matrix */
 52:   AppCtx         user;                      /* user-defined work context */
 53:   PetscInt       i,its,N,Nx = PETSC_DECIDE,Ny = PETSC_DECIDE;
 55:   PetscTruth     matrix_free = PETSC_FALSE;
 56:   PetscMPIInt    size;
 57:   PetscReal      bratu_lambda_max = 6.81,bratu_lambda_min = 0.;
 58: #if defined(PETSC_USE_LOG)
 59:   PetscLogStage  stages[2];
 60: #endif

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

 64:   PetscLogStageRegister("stage 1",&stages[0]);
 65:   PetscLogStageRegister("stage 2",&stages[1]);
 66:   for (i=0; i<2; i++) {
 67:     PetscLogStagePush(stages[i]);
 68:     user.mx = 4; user.my = 4; user.param = 6.0;
 69: 
 70:     if (i!=0) {
 71:       PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
 72:       PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
 73:       PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
 74:       if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
 75:         SETERRQ(1,"Lambda is out of range");
 76:       }
 77:     }
 78:     N = user.mx*user.my;

 80:     MPI_Comm_size(PETSC_COMM_WORLD,&size);
 81:     PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
 82:     PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
 83:     if (Nx*Ny != size && (Nx != PETSC_DECIDE || Ny != PETSC_DECIDE))
 84:       SETERRQ(1,"Incompatible number of processors:  Nx * Ny != size");
 85: 
 86:     /* Set up distributed array */
 87:     DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.mx,user.my,Nx,Ny,1,1,
 88:                       PETSC_NULL,PETSC_NULL,&user.da);
 89:     DACreateGlobalVector(user.da,&x);
 90:     VecDuplicate(x,&r);
 91:     DACreateLocalVector(user.da,&user.localX);
 92:     VecDuplicate(user.localX,&user.localF);

 94:     /* Create nonlinear solver and set function evaluation routine */
 95:     SNESCreate(PETSC_COMM_WORLD,&snes);
 96:     SNESSetType(snes,type);
 97:     SNESSetFunction(snes,r,FormFunction1,&user);

 99:     /* Set default Jacobian evaluation routine.  User can override with:
100:        -snes_mf : matrix-free Newton-Krylov method with no preconditioning
101:        (unless user explicitly sets preconditioner) 
102:        -snes_fd : default finite differencing approximation of Jacobian
103:        */
104:     PetscOptionsGetTruth(PETSC_NULL,"-snes_mf",&matrix_free,PETSC_NULL);
105:     if (!matrix_free) {
106:       PetscTruth matrix_free_operator = PETSC_FALSE;
107:       PetscOptionsGetTruth(PETSC_NULL,"-snes_mf_operator",&matrix_free_operator,PETSC_NULL);
108:       if (matrix_free_operator) matrix_free = PETSC_FALSE;
109:     }
110:     if (!matrix_free) {
111:       MatCreate(PETSC_COMM_WORLD,&J);
112:       MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,N,N);
113:       MatSetFromOptions(J);
114:       SNESSetJacobian(snes,J,J,FormJacobian1,&user);
115:     }

117:     /* Set PetscOptions, then solve nonlinear system */
118:     SNESSetFromOptions(snes);
119:     FormInitialGuess1(&user,x);
120:     SNESSolve(snes,PETSC_NULL,x);
121:     SNESGetIterationNumber(snes,&its);
122:     PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);

124:   /* Free data structures */
125:     if (!matrix_free) {
126:       MatDestroy(J);
127:     }
128:     VecDestroy(x);
129:     VecDestroy(r);
130:     VecDestroy(user.localX);
131:     VecDestroy(user.localF);
132:     SNESDestroy(snes);
133:     DADestroy(user.da);
134:   }
135:   PetscFinalize();

137:   return 0;
138: }/* --------------------  Form initial approximation ----------------- */
141: PetscErrorCode FormInitialGuess1(AppCtx *user,Vec X)
142: {
143:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,Xm,Ym,Xs,Ys;
145:   PetscReal      one = 1.0,lambda,temp1,temp,hx,hy;
146:   PetscScalar    *x;
147:   Vec            localX = user->localX;

149:   mx = user->mx;            my = user->my;            lambda = user->param;
150:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);

152:   /* Get ghost points */
153:   VecGetArray(localX,&x);
154:   temp1 = lambda/(lambda + one);
155:   DAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
156:   DAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);

158:   /* Compute initial guess */
159:   for (j=ys; j<ys+ym; j++) {
160:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
161:     for (i=xs; i<xs+xm; i++) {
162:       row = i - Xs + (j - Ys)*Xm;
163:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
164:         x[row] = 0.0;
165:         continue;
166:       }
167:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
168:     }
169:   }
170:   VecRestoreArray(localX,&x);

172:   /* Insert values into global vector */
173:   DALocalToGlobal(user->da,localX,INSERT_VALUES,X);
174:   return 0;
175: } /* --------------------  Evaluate Function F(x) --------------------- */
178: PetscErrorCode FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
179: {
180:   AppCtx         *user = (AppCtx*)ptr;
182:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym;
183:   PetscReal      two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
184:   PetscScalar    u,uxx,uyy,*x,*f;
185:   Vec            localX = user->localX,localF = user->localF;

187:   mx = user->mx;            my = user->my;            lambda = user->param;
188:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
189:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

191:   /* Get ghost points */
192:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
193:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
194:   VecGetArray(localX,&x);
195:   VecGetArray(localF,&f);
196:   DAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
197:   DAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);

199:   /* Evaluate function */
200:   for (j=ys; j<ys+ym; j++) {
201:     row = (j - Ys)*Xm + xs - Xs - 1;
202:     for (i=xs; i<xs+xm; i++) {
203:       row++;
204:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
205:         f[row] = x[row];
206:         continue;
207:       }
208:       u = x[row];
209:       uxx = (two*u - x[row-1] - x[row+1])*hydhx;
210:       uyy = (two*u - x[row-Xm] - x[row+Xm])*hxdhy;
211:       f[row] = uxx + uyy - sc*PetscExpScalar(u);
212:     }
213:   }
214:   VecRestoreArray(localX,&x);
215:   VecRestoreArray(localF,&f);

217:   /* Insert values into global vector */
218:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
219:   PetscLogFlops(11.0*ym*xm);
220:   return 0;
221: } /* --------------------  Evaluate Jacobian F'(x) --------------------- */
224: PetscErrorCode FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
225: {
226:   AppCtx         *user = (AppCtx*)ptr;
227:   Mat            jac = *J;
229:   PetscInt       i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym,col[5];
230:   PetscInt       nloc,*ltog,grow;
231:   PetscScalar    two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;
232:   Vec            localX = user->localX;

234:   mx = user->mx;            my = user->my;            lambda = user->param;
235:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
236:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

238:   /* Get ghost points */
239:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
240:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
241:   VecGetArray(localX,&x);
242:   DAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
243:   DAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);
244:   DAGetGlobalIndices(user->da,&nloc,&ltog);

246:   /* Evaluate function */
247:   for (j=ys; j<ys+ym; j++) {
248:     row = (j - Ys)*Xm + xs - Xs - 1;
249:     for (i=xs; i<xs+xm; i++) {
250:       row++;
251:       grow = ltog[row];
252:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
253:         MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
254:         continue;
255:       }
256:       v[0] = -hxdhy; col[0] = ltog[row - Xm];
257:       v[1] = -hydhx; col[1] = ltog[row - 1];
258:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
259:       v[3] = -hydhx; col[3] = ltog[row + 1];
260:       v[4] = -hxdhy; col[4] = ltog[row + Xm];
261:       MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
262:     }
263:   }
264:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
265:   VecRestoreArray(X,&x);
266:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
267:   *flag = SAME_NONZERO_PATTERN;
268:   return 0;
269: }