Actual source code: ex25.c

petsc-3.3-p7 2013-05-11
  1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods.\n";
  2: /*
  3:    u_t - alpha u_xx = A + u^2 v - (B+1) u
  4:    v_t - alpha v_xx = B u - u^2 v
  5:    0 < x < 1;
  6:    A = 1, B = 3, alpha = 1/50

  8:    Initial conditions:
  9:    u(x,0) = 1 + sin(2 pi x)
 10:    v(x,0) = 3

 12:    Boundary conditions:
 13:    u(0,t) = u(1,t) = 1
 14:    v(0,t) = v(1,t) = 3
 15: */

 17: #include <petscdmda.h>
 18: #include <petscts.h>

 20: typedef struct {
 21:   PetscScalar u,v;
 22: } Field;

 24: typedef struct _User *User;
 25: struct _User {
 26:   PetscReal A,B;                /* Reaction coefficients */
 27:   PetscReal alpha;              /* Diffusion coefficient */
 28:   PetscReal uleft,uright;       /* Dirichlet boundary conditions */
 29:   PetscReal vleft,vright;       /* Dirichlet boundary conditions */
 30: };

 32: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
 33: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 34: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
 35: static PetscErrorCode FormInitialSolution(TS,Vec,void*);

 39: int main(int argc,char **argv)
 40: {
 41:   TS                ts;         /* nonlinear solver */
 42:   Vec               X;          /* solution, residual vectors */
 43:   Mat               J;          /* Jacobian matrix */
 44:   PetscInt          steps,maxsteps,mx;
 45:   PetscErrorCode    ierr;
 46:   DM                da;
 47:   PetscReal         ftime,dt;
 48:   struct _User      user;       /* user-defined work context */
 49:   TSConvergedReason reason;

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

 53:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 54:      Create distributed array (DMDA) to manage parallel grid and vectors
 55:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 56:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,-11,2,2,PETSC_NULL,&da);

 58:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 59:      Extract global vectors from DMDA;
 60:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 61:   DMCreateGlobalVector(da,&X);

 63:   /* Initialize user application context */
 64:   PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Advection-reaction options",""); {
 65:     user.A      = 1;
 66:     user.B      = 3;
 67:     user.alpha  = 0.02;
 68:     user.uleft  = 1;
 69:     user.uright = 1;
 70:     user.vleft  = 3;
 71:     user.vright = 3;
 72:     PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,PETSC_NULL);
 73:     PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,PETSC_NULL);
 74:     PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,PETSC_NULL);
 75:     PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,PETSC_NULL);
 76:     PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,PETSC_NULL);
 77:     PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,PETSC_NULL);
 78:     PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,PETSC_NULL);
 79:   } PetscOptionsEnd();

 81:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 82:      Create timestepping solver context
 83:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 84:   TSCreate(PETSC_COMM_WORLD,&ts);
 85:   TSSetDM(ts,da);
 86:   TSSetType(ts,TSARKIMEX);
 87:   TSSetRHSFunction(ts,PETSC_NULL,FormRHSFunction,&user);
 88:   TSSetIFunction(ts,PETSC_NULL,FormIFunction,&user);
 89:   DMCreateMatrix(da,MATAIJ,&J);
 90:   TSSetIJacobian(ts,J,J,FormIJacobian,&user);

 92:   ftime = 10.0;
 93:   maxsteps = 10000;
 94:   TSSetDuration(ts,maxsteps,ftime);

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:      Set initial conditions
 98:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 99:   FormInitialSolution(ts,X,&user);
100:   TSSetSolution(ts,X);
101:   VecGetSize(X,&mx);
102:   dt   = 0.4 * user.alpha / PetscSqr(mx); /* Diffusive CFL */
103:   TSSetInitialTimeStep(ts,0.0,dt);

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:      Set runtime options
107:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108:   TSSetFromOptions(ts);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Solve nonlinear system
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   TSSolve(ts,X,&ftime);
114:   TSGetTimeStepNumber(ts,&steps);
115:   TSGetConvergedReason(ts,&reason);
116:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);

118:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119:      Free work space.
120:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121:   MatDestroy(&J);
122:   VecDestroy(&X);
123:   TSDestroy(&ts);
124:   DMDestroy(&da);
125:   PetscFinalize();
126:   return 0;
127: }

131: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
132: {
133:   User           user = (User)ptr;
134:   DM             da;
135:   DMDALocalInfo  info;
136:   PetscInt       i;
137:   Field          *x,*xdot,*f;
138:   PetscReal      hx;
139:   Vec            Xloc;

143:   TSGetDM(ts,&da);
144:   DMDAGetLocalInfo(da,&info);
145:   hx = 1.0/(PetscReal)(info.mx-1);

147:   /*
148:      Scatter ghost points to local vector,using the 2-step process
149:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
150:      By placing code between these two statements, computations can be
151:      done while messages are in transition.
152:   */
153:   DMGetLocalVector(da,&Xloc);
154:   DMGlobalToLocalBegin(da,X,INSERT_VALUES,Xloc);
155:   DMGlobalToLocalEnd(da,X,INSERT_VALUES,Xloc);

157:   /* Get pointers to vector data */
158:   DMDAVecGetArray(da,Xloc,&x);
159:   DMDAVecGetArray(da,Xdot,&xdot);
160:   DMDAVecGetArray(da,F,&f);

162:   /* Compute function over the locally owned part of the grid */
163:   for (i=info.xs; i<info.xs+info.xm; i++) {
164:     if (i == 0) {
165:       f[i].u = hx * (x[i].u - user->uleft);
166:       f[i].v = hx * (x[i].v - user->vleft);
167:     } else if (i == info.mx-1) {
168:       f[i].u = hx * (x[i].u - user->uright);
169:       f[i].v = hx * (x[i].v - user->vright);
170:     } else {
171:       f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;
172:       f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;
173:     }
174:   }

176:   /* Restore vectors */
177:   DMDAVecRestoreArray(da,Xloc,&x);
178:   DMDAVecRestoreArray(da,Xdot,&xdot);
179:   DMDAVecRestoreArray(da,F,&f);
180:   DMRestoreLocalVector(da,&Xloc);
181:   return(0);
182: }

186: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
187: {
188:   User           user = (User)ptr;
189:   DM             da;
190:   DMDALocalInfo  info;
191:   PetscInt       i;
192:   PetscReal      hx;
193:   Field          *x,*f;

197:   TSGetDM(ts,&da);
198:   DMDAGetLocalInfo(da,&info);
199:   hx = 1.0/(PetscReal)(info.mx-1);

201:   /* Get pointers to vector data */
202:   DMDAVecGetArray(da,X,&x);
203:   DMDAVecGetArray(da,F,&f);

205:   /* Compute function over the locally owned part of the grid */
206:   for (i=info.xs; i<info.xs+info.xm; i++) {
207:     PetscScalar u = x[i].u,v = x[i].v;
208:     f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
209:     f[i].v = hx*(user->B*u - u*u*v);
210:   }

212:   /* Restore vectors */
213:   DMDAVecRestoreArray(da,X,&x);
214:   DMDAVecRestoreArray(da,F,&f);
215:   return(0);
216: }

218: /* --------------------------------------------------------------------- */
219: /*
220:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
221: */
224: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ptr)
225: {
226:   User           user = (User)ptr;
228:   DMDALocalInfo  info;
229:   PetscInt       i;
230:   PetscReal      hx;
231:   DM             da;
232:   Field          *x,*xdot;

235:   TSGetDM(ts,&da);
236:   DMDAGetLocalInfo(da,&info);
237:   hx = 1.0/(PetscReal)(info.mx-1);

239:   /* Get pointers to vector data */
240:   DMDAVecGetArray(da,X,&x);
241:   DMDAVecGetArray(da,Xdot,&xdot);

243:   /* Compute function over the locally owned part of the grid */
244:   for (i=info.xs; i<info.xs+info.xm; i++) {
245:     if (i == 0 || i == info.mx-1) {
246:       const PetscInt row = i,col = i;
247:       const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
248:       MatSetValuesBlocked(*Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);
249:     } else {
250:       const PetscInt row = i,col[] = {i-1,i,i+1};
251:       const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
252:       const PetscScalar vals[2][3][2] = {{{dxxL,0},{a*hx+dxx0,0},{dxxR,0}},
253:                                          {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
254:       MatSetValuesBlocked(*Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
255:     }
256:   }

258:   /* Restore vectors */
259:   DMDAVecRestoreArray(da,X,&x);
260:   DMDAVecRestoreArray(da,Xdot,&xdot);

262:   MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);
263:   MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);
264:   if (*J != *Jpre) {
265:     MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);
266:     MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);
267:   }
268:   return(0);
269: }

273: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
274: {
275:   User user = (User)ctx;
276:   DM             da;
277:   PetscInt       i;
278:   DMDALocalInfo  info;
279:   Field          *x;
280:   PetscReal      hx;

284:   TSGetDM(ts,&da);
285:   DMDAGetLocalInfo(da,&info);
286:   hx = 1.0/(PetscReal)(info.mx-1);

288:   /* Get pointers to vector data */
289:   DMDAVecGetArray(da,X,&x);

291:   /* Compute function over the locally owned part of the grid */
292:   for (i=info.xs; i<info.xs+info.xm; i++) {
293:     PetscReal xi = i*hx;
294:     x[i].u = user->uleft*(1.-xi) + user->uright*xi + sin(2.*PETSC_PI*xi);
295:     x[i].v = user->vleft*(1.-xi) + user->vright*xi;
296:   }
297:   DMDAVecRestoreArray(da,X,&x);
298:   return(0);
299: }