Actual source code: ex25.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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 <petscdm.h>
 18: #include <petscdmda.h>
 19: #include <petscts.h>

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

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

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

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

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

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

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

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

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Create timestepping solver context
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 87:   TSCreate(PETSC_COMM_WORLD,&ts);
 88:   TSSetDM(ts,da);
 89:   TSSetType(ts,TSARKIMEX);
 90:   TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);
 91:   TSSetRHSFunction(ts,NULL,FormRHSFunction,&user);
 92:   TSSetIFunction(ts,NULL,FormIFunction,&user);
 93:   DMSetMatType(da,MATAIJ);
 94:   DMCreateMatrix(da,&J);
 95:   TSSetIJacobian(ts,J,J,FormIJacobian,&user);

 97:   ftime    = 10.0;
 98:   maxsteps = 10000;
 99:   TSSetDuration(ts,maxsteps,ftime);
100:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
101: 
102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:      Set initial conditions
104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105:   FormInitialSolution(ts,X,&user);
106:   TSSetSolution(ts,X);
107:   VecGetSize(X,&mx);
108:   hx = 1.0/(PetscReal)(mx/2-1);
109:   dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
110:   TSSetInitialTimeStep(ts,0.0,dt);

112:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113:      Set runtime options
114:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115:   TSSetFromOptions(ts);

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:      Solve nonlinear system
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120:   TSSolve(ts,X);
121:   TSGetSolveTime(ts,&ftime);
122:   TSGetTimeStepNumber(ts,&steps);
123:   TSGetConvergedReason(ts,&reason);
124:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],(double)ftime,steps);

126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:      Free work space.
128:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129:   MatDestroy(&J);
130:   VecDestroy(&X);
131:   TSDestroy(&ts);
132:   DMDestroy(&da);
133:   PetscFinalize();
134:   return 0;
135: }

139: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
140: {
141:   User           user = (User)ptr;
142:   DM             da;
143:   DMDALocalInfo  info;
144:   PetscInt       i;
145:   Field          *x,*xdot,*f;
146:   PetscReal      hx;
147:   Vec            Xloc;

151:   TSGetDM(ts,&da);
152:   DMDAGetLocalInfo(da,&info);
153:   hx   = 1.0/(PetscReal)(info.mx-1);

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

165:   /* Get pointers to vector data */
166:   DMDAVecGetArrayRead(da,Xloc,&x);
167:   DMDAVecGetArrayRead(da,Xdot,&xdot);
168:   DMDAVecGetArray(da,F,&f);

170:   /* Compute function over the locally owned part of the grid */
171:   for (i=info.xs; i<info.xs+info.xm; i++) {
172:     if (i == 0) {
173:       f[i].u = hx * (x[i].u - user->uleft);
174:       f[i].v = hx * (x[i].v - user->vleft);
175:     } else if (i == info.mx-1) {
176:       f[i].u = hx * (x[i].u - user->uright);
177:       f[i].v = hx * (x[i].v - user->vright);
178:     } else {
179:       f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;
180:       f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;
181:     }
182:   }

184:   /* Restore vectors */
185:   DMDAVecRestoreArrayRead(da,Xloc,&x);
186:   DMDAVecRestoreArrayRead(da,Xdot,&xdot);
187:   DMDAVecRestoreArray(da,F,&f);
188:   DMRestoreLocalVector(da,&Xloc);
189:   return(0);
190: }

194: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
195: {
196:   User           user = (User)ptr;
197:   DM             da;
198:   DMDALocalInfo  info;
199:   PetscInt       i;
200:   PetscReal      hx;
201:   Field          *x,*f;

205:   TSGetDM(ts,&da);
206:   DMDAGetLocalInfo(da,&info);
207:   hx   = 1.0/(PetscReal)(info.mx-1);

209:   /* Get pointers to vector data */
210:   DMDAVecGetArrayRead(da,X,&x);
211:   DMDAVecGetArray(da,F,&f);

213:   /* Compute function over the locally owned part of the grid */
214:   for (i=info.xs; i<info.xs+info.xm; i++) {
215:     PetscScalar u = x[i].u,v = x[i].v;
216:     f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
217:     f[i].v = hx*(user->B*u - u*u*v);
218:   }

220:   /* Restore vectors */
221:   DMDAVecRestoreArrayRead(da,X,&x);
222:   DMDAVecRestoreArray(da,F,&f);
223:   return(0);
224: }

226: /* --------------------------------------------------------------------- */
227: /*
228:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
229: */
232: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
233: {
234:   User           user = (User)ptr;
236:   DMDALocalInfo  info;
237:   PetscInt       i;
238:   PetscReal      hx;
239:   DM             da;
240:   Field          *x,*xdot;

243:   TSGetDM(ts,&da);
244:   DMDAGetLocalInfo(da,&info);
245:   hx   = 1.0/(PetscReal)(info.mx-1);

247:   /* Get pointers to vector data */
248:   DMDAVecGetArrayRead(da,X,&x);
249:   DMDAVecGetArrayRead(da,Xdot,&xdot);

251:   /* Compute function over the locally owned part of the grid */
252:   for (i=info.xs; i<info.xs+info.xm; i++) {
253:     if (i == 0 || i == info.mx-1) {
254:       const PetscInt    row        = i,col = i;
255:       const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
256:       MatSetValuesBlocked(Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);
257:     } else {
258:       const PetscInt    row           = i,col[] = {i-1,i,i+1};
259:       const PetscScalar dxxL          = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
260:       const PetscScalar vals[2][3][2] = {{{dxxL,0},{a *hx+dxx0,0},{dxxR,0}},
261:                                          {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
262:       MatSetValuesBlocked(Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);
263:     }
264:   }

266:   /* Restore vectors */
267:   DMDAVecRestoreArrayRead(da,X,&x);
268:   DMDAVecRestoreArrayRead(da,Xdot,&xdot);

270:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
271:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
272:   if (J != Jpre) {
273:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
274:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
275:   }
276:   return(0);
277: }

281: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
282: {
283:   User           user = (User)ctx;
284:   DM             da;
285:   PetscInt       i;
286:   DMDALocalInfo  info;
287:   Field          *x;
288:   PetscReal      hx;

292:   TSGetDM(ts,&da);
293:   DMDAGetLocalInfo(da,&info);
294:   hx   = 1.0/(PetscReal)(info.mx-1);

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

299:   /* Compute function over the locally owned part of the grid */
300:   for (i=info.xs; i<info.xs+info.xm; i++) {
301:     PetscReal xi = i*hx;
302:     x[i].u = user->uleft*(1.-xi) + user->uright*xi + PetscSinReal(2.*PETSC_PI*xi);
303:     x[i].v = user->vleft*(1.-xi) + user->vright*xi;
304:   }
305:   DMDAVecRestoreArray(da,X,&x);
306:   return(0);
307: }