Actual source code: ex9.c

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

  3: static char help[] = "This program demonstrates use of the SNES package. Solve systems of\n\
  4: nonlinear equations in parallel.  This example uses matrix-free Krylov\n\
  5: Newton methods with no preconditioner to solve the Bratu (SFI - solid fuel\n\
  6: ignition) test problem. The command line options are:\n\
  7:    -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\
  8:       problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)\n\
  9:    -mx <xg>, where <xg> = number of grid points in the x-direction\n\
 10:    -my <yg>, where <yg> = number of grid points in the y-direction\n\
 11:    -mz <zg>, where <zg> = number of grid points in the z-direction\n\n";

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

 28:  #include petscsnes.h
 29:  #include petscda.h

 31: typedef struct {
 32:     PetscReal param;           /* test problem nonlinearity parameter */
 33:     int       mx,my,mz;      /* discretization in x,y,z-directions */
 34:     Vec       localX,localF;  /* ghosted local vectors */
 35:     DA        da;              /* distributed array datastructure */
 36: } AppCtx;

 38: extern int FormFunction1(SNES,Vec,Vec,void*),FormInitialGuess1(AppCtx*,Vec);

 42: int main(int argc,char **argv)
 43: {
 44:   SNES          snes;                 /* nonlinear solver */
 45:   KSP           ksp;                 /* linear solver */
 46:   PC            pc;                   /* preconditioner */
 47:   Mat           J;                    /* Jacobian matrix */
 48:   AppCtx        user;                 /* user-defined application context */
 49:   Vec           x,r;                  /* vectors */
 50:   DAStencilType stencil = DA_STENCIL_BOX;
 51:   int           ierr,its;
 52:   PetscTruth    flg;
 53:   int           Nx = PETSC_DECIDE,Ny = PETSC_DECIDE,Nz = PETSC_DECIDE;
 54:   PetscReal     bratu_lambda_max = 6.81,bratu_lambda_min = 0.;

 56:   PetscInitialize(&argc,&argv,(char *)0,help);
 57:   PetscOptionsHasName(PETSC_NULL,"-star",&flg);
 58:   if (flg) stencil = DA_STENCIL_STAR;

 60:   user.mx    = 4;
 61:   user.my    = 4;
 62:   user.mz    = 4;
 63:   user.param = 6.0;
 64:   PetscOptionsGetInt(PETSC_NULL,"-mx",&user.mx,PETSC_NULL);
 65:   PetscOptionsGetInt(PETSC_NULL,"-my",&user.my,PETSC_NULL);
 66:   PetscOptionsGetInt(PETSC_NULL,"-mz",&user.mz,PETSC_NULL);
 67:   PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
 68:   PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
 69:   PetscOptionsGetInt(PETSC_NULL,"-Nz",&Nz,PETSC_NULL);
 70:   PetscOptionsGetReal(PETSC_NULL,"-par",&user.param,PETSC_NULL);
 71:   if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) {
 72:     SETERRQ(1,"Lambda is out of range");
 73:   }
 74: 
 75:   /* Set up distributed array */
 76:   DACreate3d(PETSC_COMM_WORLD,DA_NONPERIODIC,stencil,user.mx,user.my,user.mz,
 77:                     Nx,Ny,Nz,1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&user.da);
 78:   DACreateGlobalVector(user.da,&x);
 79:   VecDuplicate(x,&r);
 80:   DACreateLocalVector(user.da,&user.localX);
 81:   VecDuplicate(user.localX,&user.localF);

 83:   /* Create nonlinear solver */
 84:   SNESCreate(PETSC_COMM_WORLD,&snes);
 85:   /* Set various routines and options */
 86:   SNESSetFunction(snes,r,FormFunction1,(void *)&user);
 87:   MatCreateSNESMF(snes,x,&J);
 88:   SNESSetJacobian(snes,J,J,MatSNESMFComputeJacobian,&user);
 89:   SNESSetFromOptions(snes);

 91:   /* Force no preconditioning to be used.  Note that this overrides whatever
 92:      choices may have been specified in the options database. */
 93:   SNESGetKSP(snes,&ksp);
 94:   KSPGetPC(ksp,&pc);
 95:   PCSetType(pc,PCNONE);

 97:   /* Solve nonlinear system */
 98:   FormInitialGuess1(&user,x);
 99:   SNESSolve(snes,x);
100:   SNESGetIterationNumber(snes,&its);
101:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %d\n",its);

103:   /* Free data structures */
104:   VecDestroy(user.localX);
105:   VecDestroy(user.localF);
106:   DADestroy(user.da);
107:   VecDestroy(x); VecDestroy(r);
108:   MatDestroy(J); SNESDestroy(snes);

110:   PetscFinalize();
111:   return 0;
112: }/* --------------------  Form initial approximation ----------------- */
115: int FormInitialGuess1(AppCtx *user,Vec X)
116: {
117:   int          i,j,k,loc,mx,my,mz,ierr,xs,ys,zs,xm,ym,zm,Xm,Ym,Zm,Xs,Ys,Zs,base1;
118:   PetscReal    one = 1.0,lambda,temp1,temp,Hx,Hy;
119:   PetscScalar  *x;
120:   Vec          localX = user->localX;

122:   mx         = user->mx; my         = user->my; mz = user->mz; lambda = user->param;
123:   Hx     = one / (PetscReal)(mx-1);     Hy     = one / (PetscReal)(my-1);

125:   VecGetArray(localX,&x);
126:   temp1 = lambda/(lambda + one);
127:   DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
128:   DAGetGhostCorners(user->da,&Xs,&Ys,&Zs,&Xm,&Ym,&Zm);
129: 
130:   for (k=zs; k<zs+zm; k++) {
131:     base1 = (Xm*Ym)*(k-Zs);
132:     for (j=ys; j<ys+ym; j++) {
133:       temp = (PetscReal)(PetscMin(j,my-j-1))*Hy;
134:       for (i=xs; i<xs+xm; i++) {
135:         loc = base1 + i-Xs + (j-Ys)*Xm;
136:         if (i == 0 || j == 0 || k == 0 || i==mx-1 || j==my-1 || k==mz-1) {
137:           x[loc] = 0.0;
138:           continue;
139:         }
140:         x[loc] = temp1*sqrt(PetscMin((PetscReal)(PetscMin(i,mx-i-1))*Hx,temp));
141:       }
142:     }
143:   }

145:   VecRestoreArray(localX,&x);
146:   /* stick values into global vector */
147:   DALocalToGlobal(user->da,localX,INSERT_VALUES,X);
148:   return 0;
149: }/* --------------------  Evaluate Function F(x) --------------------- */
152: int FormFunction1(SNES snes,Vec X,Vec F,void *ptr)
153: {
154:   AppCtx       *user = (AppCtx*)ptr;
155:   int          ierr,i,j,k,loc,mx,my,mz,xs,ys,zs,xm,ym,zm,Xs,Ys,Zs,Xm,Ym,Zm;
156:   int          base1,base2;
157:   PetscReal    two = 2.0,one = 1.0,lambda,Hx,Hy,Hz,HxHzdHy,HyHzdHx,HxHydHz;
158:   PetscScalar  u,uxx,uyy,sc,*x,*f,uzz;
159:   Vec          localX = user->localX,localF = user->localF;

161:   mx      = user->mx; my = user->my; mz = user->mz; lambda = user->param;
162:   Hx      = one / (PetscReal)(mx-1);
163:   Hy      = one / (PetscReal)(my-1);
164:   Hz      = one / (PetscReal)(mz-1);
165:   sc      = Hx*Hy*Hz*lambda; HxHzdHy  = Hx*Hz/Hy; HyHzdHx  = Hy*Hz/Hx;
166:   HxHydHz = Hx*Hy/Hz;

168:   DAGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
169:   DAGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
170:   VecGetArray(localX,&x);
171:   VecGetArray(localF,&f);

173:   DAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
174:   DAGetGhostCorners(user->da,&Xs,&Ys,&Zs,&Xm,&Ym,&Zm);

176:   for (k=zs; k<zs+zm; k++) {
177:     base1 = (Xm*Ym)*(k-Zs);
178:     for (j=ys; j<ys+ym; j++) {
179:       base2 = base1 + (j-Ys)*Xm;
180:       for (i=xs; i<xs+xm; i++) {
181:         loc = base2 + (i-Xs);
182:         if (i == 0 || j == 0 || k== 0 || i == mx-1 || j == my-1 || k == mz-1) {
183:           f[loc] = x[loc];
184:         }
185:         else {
186:           u = x[loc];
187:           uxx = (two*u - x[loc-1] - x[loc+1])*HyHzdHx;
188:           uyy = (two*u - x[loc-Xm] - x[loc+Xm])*HxHzdHy;
189:           uzz = (two*u - x[loc-Xm*Ym] - x[loc+Xm*Ym])*HxHydHz;
190:           f[loc] = uxx + uyy + uzz - sc*PetscExpScalar(u);
191:         }
192:       }
193:     }
194:   }
195:   VecRestoreArray(localX,&x);
196:   VecRestoreArray(localF,&f);
197:   /* stick values into global vector */
198:   DALocalToGlobal(user->da,localF,INSERT_VALUES,F);
199:   PetscLogFlops(11*ym*xm*zm);
200:   return 0;
201: }
202: 




207: