Actual source code: ex9.c

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

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

 27: #include <petscsnes.h>
 28: #include <petscdmda.h>

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


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

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

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

 80:   /* Create nonlinear solver */
 81:   SNESCreate(PETSC_COMM_WORLD,&snes);
 82:   /* Set various routines and options */
 83:   SNESSetFunction(snes,r,FormFunction1,(void*)&user);
 84:   MatCreateSNESMF(snes,&J);
 85:   SNESSetJacobian(snes,J,J,MatMFFDComputeJacobian,&user);
 86:   SNESSetFromOptions(snes);

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

 94:   /* Solve nonlinear system */
 95:   FormInitialGuess1(&user,x);
 96:   SNESSolve(snes,PETSC_NULL,x);
 97:   SNESGetIterationNumber(snes,&its);
 98:   PetscPrintf(PETSC_COMM_WORLD,"Number of Newton iterations = %D\n",its);

100:   /* Free data structures */
101:   VecDestroy(&user.localX);
102:   VecDestroy(&user.localF);
103:   DMDestroy(&user.da);
104:   VecDestroy(&x); VecDestroy(&r);
105:   MatDestroy(&J); SNESDestroy(&snes);

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

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

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

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

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

167:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
168:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
169:   VecGetArray(localX,&x);
170:   VecGetArray(localF,&f);

172:   DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);
173:   DMDAGetGhostCorners(user->da,&Xs,&Ys,&Zs,&Xm,&Ym,&Zm);

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




207: