Actual source code: ex25.c

  1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d formulated as a PDAE. Demonstrates solving PDEs with algebraic constraints (PDAE).\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 *);

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

 50:   PetscInitialize(&argc, &argv, (char *)0, help);
 51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 52:      Create distributed array (DMDA) to manage parallel grid and vectors
 53:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 54:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da);
 55:   DMSetFromOptions(da);
 56:   DMSetUp(da);

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

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

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

 96:   ftime = 10.0;
 97:   TSSetMaxTime(ts, ftime);
 98:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Set initial conditions
102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   FormInitialSolution(ts, X, &user);
104:   TSSetSolution(ts, X);
105:   VecGetSize(X, &mx);
106:   hx = 1.0 / (PetscReal)(mx / 2 - 1);
107:   dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
108:   TSSetTimeStep(ts, dt);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Set runtime options
112:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   TSSetFromOptions(ts);

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

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

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

146:   TSGetDM(ts, &da);
147:   DMDAGetLocalInfo(da, &info);
148:   hx = 1.0 / (PetscReal)(info.mx - 1);

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

160:   /* Get pointers to vector data */
161:   DMDAVecGetArrayRead(da, Xloc, &x);
162:   DMDAVecGetArrayRead(da, Xdot, &xdot);
163:   DMDAVecGetArray(da, F, &f);

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

179:   /* Restore vectors */
180:   DMDAVecRestoreArrayRead(da, Xloc, &x);
181:   DMDAVecRestoreArrayRead(da, Xdot, &xdot);
182:   DMDAVecRestoreArray(da, F, &f);
183:   DMRestoreLocalVector(da, &Xloc);
184:   return 0;
185: }

187: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
188: {
189:   User          user = (User)ptr;
190:   DM            da;
191:   DMDALocalInfo info;
192:   PetscInt      i;
193:   PetscReal     hx;
194:   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:   DMDAVecGetArrayRead(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:   DMDAVecRestoreArrayRead(da, X, &x);
214:   DMDAVecRestoreArray(da, F, &f);
215:   return 0;
216: }

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

232:   TSGetDM(ts, &da);
233:   DMDAGetLocalInfo(da, &info);
234:   hx = 1.0 / (PetscReal)(info.mx - 1);

236:   /* Get pointers to vector data */
237:   DMDAVecGetArrayRead(da, X, &x);
238:   DMDAVecGetArrayRead(da, Xdot, &xdot);

240:   /* Compute function over the locally owned part of the grid */
241:   for (i = info.xs; i < info.xs + info.xm; i++) {
242:     if (i == 0 || i == info.mx - 1) {
243:       const PetscInt    row = i, col = i;
244:       const PetscScalar vals[2][2] = {
245:         {hx, 0 },
246:         {0,  hx}
247:       };
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] = {
253:         {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
254:         {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
255:       };
256:       MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
257:     }
258:   }

260:   /* Restore vectors */
261:   DMDAVecRestoreArrayRead(da, X, &x);
262:   DMDAVecRestoreArrayRead(da, Xdot, &xdot);

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

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;

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

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

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

300: /*TEST

302:     test:
303:       args: -nox -da_grid_x 20 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none

305: TEST*/