Actual source code: ex7.c

  1: /*$Id: ex7.c,v 1.7 2001/08/21 21:04:28 bsmith Exp $*/

  3: /* Program usage:  mpirun -np <procs> ex5 [-help] [all PETSc options] */

  5: static char help[] = "Nonlinear, time-dependent PDE in 2d.n";


  8: /* 
  9:    Include "petscda.h" so that we can use distributed arrays (DAs).
 10:    Include "petscts.h" so that we can use SNES solvers.  Note that this
 11:    file automatically includes:
 12:      petsc.h       - base PETSc routines   petscvec.h - vectors
 13:      petscsys.h    - system routines       petscmat.h - matrices
 14:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 15:      petscviewer.h - viewers               petscpc.h  - preconditioners
 16:      petscsles.h   - linear solvers
 17: */
 18:  #include petscda.h
 19:  #include petscts.h


 22: /* 
 23:    User-defined routines
 24: */
 25: extern int FormFunction(TS,PetscReal,Vec,Vec,void*),FormInitialSolution(DA,Vec);
 26: extern int Monitor(TS,int,PetscReal,Vec,void*);

 28: #undef __FUNCT__
 30: int main(int argc,char **argv)
 31: {
 32:   TS                     ts;                 /* nonlinear solver */
 33:   Vec                    x,r;                  /* solution, residual vectors */
 34:   Mat                    J;                    /* Jacobian matrix */
 35:   int                    steps;                  /* iterations for convergence */
 36:   int                    ierr;
 37:   DA                     da;
 38:   MatFDColoring          matfdcoloring;
 39:   ISColoring             iscoloring;
 40:   PetscReal              ftime;

 42:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 43:      Initialize program
 44:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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


 49:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 50:      Create distributed array (DA) to manage parallel grid and vectors
 51:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 52:   DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,8,8,PETSC_DECIDE,PETSC_DECIDE,
 53:                     1,1,PETSC_NULL,PETSC_NULL,&da);

 55:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 56:      Extract global vectors from DA; then duplicate for remaining
 57:      vectors that are the same types
 58:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 59:   DACreateGlobalVector(da,&x);
 60:   VecDuplicate(x,&r);

 62:   TSCreate(PETSC_COMM_WORLD,&ts);
 63:   TSSetProblemType(ts,TS_NONLINEAR);
 64:   TSSetRHSFunction(ts,FormFunction,da);


 67:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 68:      Create matrix data structure; set Jacobian evaluation routine

 70:      Set Jacobian matrix data structure and default Jacobian evaluation
 71:      routine. User can override with:
 72:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
 73:                 (unless user explicitly sets preconditioner) 
 74:      -snes_mf_operator : form preconditioning matrix as set by the user,
 75:                          but use matrix-free approx for Jacobian-vector
 76:                          products within Newton-Krylov method

 78:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 79:   DAGetColoring(da,IS_COLORING_LOCAL,&iscoloring);
 80:   DAGetMatrix(da,MATMPIAIJ,&J);
 81:   MatFDColoringCreate(J,iscoloring,&matfdcoloring);
 82:   ISColoringDestroy(iscoloring);
 83:   MatFDColoringSetFunction(matfdcoloring,(int (*)(void))FormFunction,da);
 84:   MatFDColoringSetFromOptions(matfdcoloring);
 85:   TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
 86:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87:      Customize nonlinear solver; set runtime options
 88:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 89:   TSSetInitialTimeStep(ts,0.0,.0001);
 90:   TSSetType(ts,TS_BEULER);
 91:   TSSetDuration(ts,100,1.0);
 92:   TSSetFromOptions(ts);
 93:   TSSetMonitor(ts,Monitor,0,0);

 95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 96:      Set initial conditions
 97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 98:   FormInitialSolution(da,x);

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Solve nonlinear system
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   TSSetSolution(ts,x);
104:   TSStep(ts,&steps,&ftime);


107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Free work space.  All PETSc objects should be destroyed when they
109:      are no longer needed.
110:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

112:   MatDestroy(J);
113:   MatFDColoringDestroy(matfdcoloring);
114:   VecDestroy(x);
115:   VecDestroy(r);
116:   TSDestroy(ts);
117:   DADestroy(da);
118:   PetscFinalize();

120:   return(0);
121: }
122: /* ------------------------------------------------------------------- */
123: #undef __FUNCT__
125: /* 
126:    FormFunction - Evaluates nonlinear function, F(x).

128:    Input Parameters:
129: .  ts - the TS context
130: .  X - input vector
131: .  ptr - optional user-defined context, as set by SNESSetFunction()

133:    Output Parameter:
134: .  F - function vector
135:  */
136: int FormFunction(TS ts,PetscReal time,Vec X,Vec F,void *ptr)
137: {
138:   DA           da = (DA)ptr;
139:   int          ierr,i,j,Mx,My,xs,ys,xm,ym;
140:   PetscReal    two = 2.0,hx,hy,hxdhy,hydhx,sx,sy;
141:   PetscScalar  u,uxx,uyy,**x,**f;
142:   Vec          localX;

145:   DAGetLocalVector(da,&localX);
146:   DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
147:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

149:   hx     = 1.0/(PetscReal)(Mx-1); sx = 1.0/(hx*hx);
150:   hy     = 1.0/(PetscReal)(My-1); sy = 1.0/(hy*hy);
151:   hxdhy  = hx/hy;
152:   hydhx  = hy/hx;

154:   /*
155:      Scatter ghost points to local vector,using the 2-step process
156:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
157:      By placing code between these two statements, computations can be
158:      done while messages are in transition.
159:   */
160:   DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX);
161:   DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX);

163:   /*
164:      Get pointers to vector data
165:   */
166:   DAVecGetArray(da,localX,(void**)&x);
167:   DAVecGetArray(da,F,(void**)&f);

169:   /*
170:      Get local grid boundaries
171:   */
172:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

174:   /*
175:      Compute function over the locally owned part of the grid
176:   */
177:   for (j=ys; j<ys+ym; j++) {
178:     for (i=xs; i<xs+xm; i++) {
179:       if (i == 0 || j == 0 || i == Mx-1 || j == My-1) {
180:         f[j][i] = x[j][i];
181:         continue;
182:       }
183:       u       = x[j][i];
184:       uxx     = (two*u - x[j][i-1] - x[j][i+1])*sx;
185:       uyy     = (two*u - x[j-1][i] - x[j+1][i])*sy;
186:       /*      f[j][i] = -(uxx + uyy); */
187:       f[j][i] = -u*(uxx + uyy) - (4.0 - 1.0)*((x[j][i+1] - x[j][i-1])*(x[j][i+1] - x[j][i-1])*.25*sx +
188:                                             (x[j+1][i] - x[j-1][i])*(x[j+1][i] - x[j-1][i])*.25*sy);
189:     }
190:   }

192:   /*
193:      Restore vectors
194:   */
195:   DAVecRestoreArray(da,localX,(void**)&x);
196:   DAVecRestoreArray(da,F,(void**)&f);
197:   DARestoreLocalVector(da,&localX);
198:   PetscLogFlops(11*ym*xm);
199:   return(0);
200: }

202: /* ------------------------------------------------------------------- */
203: #undef __FUNCT__
205: int FormInitialSolution(DA da,Vec U)
206: {
207:   int         ierr,i,j,xs,ys,xm,ym,Mx,My;
208:   PetscScalar **u;
209:   PetscReal   hx,hy,x,y,r;

212:   DAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
213:                    PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);

215:   hx     = 1.0/(PetscReal)(Mx-1);
216:   hy     = 1.0/(PetscReal)(My-1);

218:   /*
219:      Get pointers to vector data
220:   */
221:   DAVecGetArray(da,U,(void**)&u);

223:   /*
224:      Get local grid boundaries
225:   */
226:   DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);

228:   /*
229:      Compute function over the locally owned part of the grid
230:   */
231:   for (j=ys; j<ys+ym; j++) {
232:     y = j*hy;
233:     for (i=xs; i<xs+xm; i++) {
234:       x = i*hx;
235:       r = PetscSqrtScalar((x-.5)*(x-.5) + (y-.5)*(y-.5));
236:       if (r < .125) {
237:         u[j][i] = PetscExpScalar(-30.0*r*r*r);
238:       } else {
239:         u[j][i] = 0.0;
240:       }
241:     }
242:   }

244:   /*
245:      Restore vectors
246:   */
247:   DAVecRestoreArray(da,U,(void**)&u);
248:   return(0);
249: }

251: #undef __FUNCT__  
253: int Monitor(TS ts,int step,PetscReal ptime,Vec v,void *ctx)
254: {
255:   int       ierr;
256:   PetscReal norm;
257:   MPI_Comm  comm;

260:   VecNorm(v,NORM_2,&norm);
261:   PetscObjectGetComm((PetscObject)ts,&comm);
262:   PetscPrintf(comm,"timestep %d time %g norm %gn",step,ptime,norm);
263:   return(0);
264: }