Actual source code: ex13.c

  1: /*$Id: ex13.c,v 1.33 2001/08/07 21:31:12 bsmith Exp $*/

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

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

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

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

 42: extern int FormFunction1(SNES,Vec,Vec,void*),FormInitialGuess1(AppCtx*,Vec);
 43: extern int FormJacobian1(SNES,Vec,Mat*,Mat*,MatStructure*,void*);

 47: int main(int argc,char **argv)
 48: {
 49:   SNES          snes;                      /* nonlinear solver */
 50:   SNESType      type = SNESLS;             /* nonlinear solution method */
 51:   Vec           x,r;                       /* solution, residual vectors */
 52:   Mat           J;                         /* Jacobian matrix */
 53:   AppCtx        user;                      /* user-defined work context */
 54:   int           i,ierr,its,N,Nx = PETSC_DECIDE,Ny = PETSC_DECIDE;
 55:   PetscTruth    matrix_free;
 56:   int           size;
 57:   PetscReal     bratu_lambda_max = 6.81,bratu_lambda_min = 0.;

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

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

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

 89:     /* Create nonlinear solver and set function evaluation routine */
 90:     SNESCreate(PETSC_COMM_WORLD,&snes);
 91:     SNESSetType(snes,type);
 92:     SNESSetFunction(snes,r,FormFunction1,&user);

 94:     /* Set default Jacobian evaluation routine.  User can override with:
 95:        -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 96:        (unless user explicitly sets preconditioner) 
 97:        -snes_fd : default finite differencing approximation of Jacobian
 98:        */
 99:     PetscOptionsHasName(PETSC_NULL,"-snes_mf",&matrix_free);
100:     if (!matrix_free) {
101:       MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,N,N,&J);
102:       MatSetFromOptions(J);
103:       SNESSetJacobian(snes,J,J,FormJacobian1,&user);
104:     }

106:     /* Set PetscOptions, then solve nonlinear system */
107:     SNESSetFromOptions(snes);
108:     FormInitialGuess1(&user,x);
109:     SNESSolve(snes,x,&its);
110:     PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n",its);

112:   /* Free data structures */
113:     if (!matrix_free) {
114:       MatDestroy(J);
115:     }
116:     VecDestroy(x);
117:     VecDestroy(r);
118:     VecDestroy(user.localX);
119:     VecDestroy(user.localF);
120:     SNESDestroy(snes);
121:     DADestroy(user.da);
122:   }
123:   PetscFinalize();

125:   return 0;
126: }/* --------------------  Form initial approximation ----------------- */
129: int FormInitialGuess1(AppCtx *user,Vec X)
130: {
131:   int          i,j,row,mx,my,ierr,xs,ys,xm,ym,Xm,Ym,Xs,Ys;
132:   PetscReal    one = 1.0,lambda,temp1,temp,hx,hy;
133:   PetscScalar  *x;
134:   Vec          localX = user->localX;

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

139:   /* Get ghost points */
140:   VecGetArray(localX,&x);
141:   temp1 = lambda/(lambda + one);
142:   DAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
143:   DAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);

145:   /* Compute initial guess */
146:   for (j=ys; j<ys+ym; j++) {
147:     temp = (PetscReal)(PetscMin(j,my-j-1))*hy;
148:     for (i=xs; i<xs+xm; i++) {
149:       row = i - Xs + (j - Ys)*Xm;
150:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
151:         x[row] = 0.0;
152:         continue;
153:       }
154:       x[row] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*hx,temp));
155:     }
156:   }
157:   VecRestoreArray(localX,&x);

159:   /* Insert values into global vector */
160:   DALocalToGlobal(user->da,localX,INSERT_VALUES,X);
161:   return 0;
162: } /* --------------------  Evaluate Function F(x) --------------------- */
165: int FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
166: {
167:   AppCtx       *user = (AppCtx*)ptr;
168:   int          ierr,i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym;
169:   PetscReal    two = 2.0,one = 1.0,lambda,hx,hy,hxdhy,hydhx,sc;
170:   PetscScalar  u,uxx,uyy,*x,*f;
171:   Vec          localX = user->localX,localF = user->localF;

173:   mx = user->mx;            my = user->my;            lambda = user->param;
174:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
175:   sc = hx*hy*lambda;        hxdhy = hx/hy;            hydhx = hy/hx;

177:   /* Get ghost points */
178:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
179:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
180:   VecGetArray(localX,&x);
181:   VecGetArray(localF,&f);
182:   DAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
183:   DAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);

185:   /* Evaluate function */
186:   for (j=ys; j<ys+ym; j++) {
187:     row = (j - Ys)*Xm + xs - Xs - 1;
188:     for (i=xs; i<xs+xm; i++) {
189:       row++;
190:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
191:         f[row] = x[row];
192:         continue;
193:       }
194:       u = x[row];
195:       uxx = (two*u - x[row-1] - x[row+1])*hydhx;
196:       uyy = (two*u - x[row-Xm] - x[row+Xm])*hxdhy;
197:       f[row] = uxx + uyy - sc*PetscExpScalar(u);
198:     }
199:   }
200:   VecRestoreArray(localX,&x);
201:   VecRestoreArray(localF,&f);

203:   /* Insert values into global vector */
204:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
205:   PetscLogFlops(11*ym*xm);
206:   return 0;
207: } /* --------------------  Evaluate Jacobian F'(x) --------------------- */
210: int FormJacobian1(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
211: {
212:   AppCtx  *user = (AppCtx*)ptr;
213:   Mat     jac = *J;
214:   int     ierr,i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym,col[5];
215:   int     nloc,*ltog,grow;
216:   PetscScalar  two = 2.0,one = 1.0,lambda,v[5],hx,hy,hxdhy,hydhx,sc,*x;
217:   Vec     localX = user->localX;

219:   mx = user->mx;            my = user->my;            lambda = user->param;
220:   hx = one/(PetscReal)(mx-1);  hy = one/(PetscReal)(my-1);
221:   sc = hx*hy;               hxdhy = hx/hy;            hydhx = hy/hx;

223:   /* Get ghost points */
224:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
225:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
226:   VecGetArray(localX,&x);
227:   DAGetCorners(user->da,&xs,&ys,0,&xm,&ym,0);
228:   DAGetGhostCorners(user->da,&Xs,&Ys,0,&Xm,&Ym,0);
229:   DAGetGlobalIndices(user->da,&nloc,&ltog);

231:   /* Evaluate function */
232:   for (j=ys; j<ys+ym; j++) {
233:     row = (j - Ys)*Xm + xs - Xs - 1;
234:     for (i=xs; i<xs+xm; i++) {
235:       row++;
236:       grow = ltog[row];
237:       if (i == 0 || j == 0 || i == mx-1 || j == my-1) {
238:         MatSetValues(jac,1,&grow,1,&grow,&one,INSERT_VALUES);
239:         continue;
240:       }
241:       v[0] = -hxdhy; col[0] = ltog[row - Xm];
242:       v[1] = -hydhx; col[1] = ltog[row - 1];
243:       v[2] = two*(hydhx + hxdhy) - sc*lambda*PetscExpScalar(x[row]); col[2] = grow;
244:       v[3] = -hydhx; col[3] = ltog[row + 1];
245:       v[4] = -hxdhy; col[4] = ltog[row + Xm];
246:       MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
247:     }
248:   }
249:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
250:   VecRestoreArray(X,&x);
251:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
252:   *flag = SAME_NONZERO_PATTERN;
253:   return 0;
254: }