Actual source code: ex7.c


  2: static char help[] = "Time-dependent PDE in 2d for calculating joint PDF. \n";
  3: /*
  4:    p_t = -x_t*p_x -y_t*p_y + f(t)*p_yy
  5:    xmin < x < xmax, ymin < y < ymax;

  7:    Boundary conditions Neumman using mirror values

  9:    Note that x_t and y_t in the above are given functions of x and y; they are not derivatives of x and y.
 10:    x_t = (y - ws)  y_t = (ws/2H)*(Pm - Pmax*sin(x))

 12: */

 14: #include <petscdm.h>
 15: #include <petscdmda.h>
 16: #include <petscts.h>

 18: /*
 19:    User-defined data structures and routines
 20: */
 21: typedef struct {
 22:   PetscScalar    ws;         /* Synchronous speed */
 23:   PetscScalar    H;          /* Inertia constant */
 24:   PetscScalar    D;          /* Damping constant */
 25:   PetscScalar    Pmax;       /* Maximum power output of generator */
 26:   PetscScalar    PM_min;     /* Mean mechanical power input */
 27:   PetscScalar    lambda;     /* correlation time */
 28:   PetscScalar    q;          /* noise strength */
 29:   PetscScalar    mux;        /* Initial average angle */
 30:   PetscScalar    sigmax;     /* Standard deviation of initial angle */
 31:   PetscScalar    muy;        /* Average speed */
 32:   PetscScalar    sigmay;     /* standard deviation of initial speed */
 33:   PetscScalar    rho;        /* Cross-correlation coefficient */
 34:   PetscScalar    xmin;       /* left boundary of angle */
 35:   PetscScalar    xmax;       /* right boundary of angle */
 36:   PetscScalar    ymin;       /* bottom boundary of speed */
 37:   PetscScalar    ymax;       /* top boundary of speed */
 38:   PetscScalar    dx;         /* x step size */
 39:   PetscScalar    dy;         /* y step size */
 40:   PetscScalar    disper_coe; /* Dispersion coefficient */
 41:   DM             da;
 42:   PetscInt       st_width; /* Stencil width */
 43:   DMBoundaryType bx;       /* x boundary type */
 44:   DMBoundaryType by;       /* y boundary type */
 45:   PetscBool      nonoiseinitial;
 46: } AppCtx;

 48: PetscErrorCode Parameter_settings(AppCtx *);
 49: PetscErrorCode ini_bou(Vec, AppCtx *);
 50: PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
 51: PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
 52: PetscErrorCode PostStep(TS);

 54: int main(int argc, char **argv)
 55: {
 56:   Vec         x;    /* Solution vector */
 57:   TS          ts;   /* Time-stepping context */
 58:   AppCtx      user; /* Application context */
 59:   PetscViewer viewer;

 62:   PetscInitialize(&argc, &argv, "petscopt_ex7", help);

 64:   /* Get physics and time parameters */
 65:   Parameter_settings(&user);
 66:   /* Create a 2D DA with dof = 1 */
 67:   DMDACreate2d(PETSC_COMM_WORLD, user.bx, user.by, DMDA_STENCIL_STAR, 4, 4, PETSC_DECIDE, PETSC_DECIDE, 1, user.st_width, NULL, NULL, &user.da);
 68:   DMSetFromOptions(user.da);
 69:   DMSetUp(user.da);
 70:   /* Set x and y coordinates */
 71:   DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0, 0);
 72:   DMDASetCoordinateName(user.da, 0, "X - the angle");
 73:   DMDASetCoordinateName(user.da, 1, "Y - the speed");

 75:   /* Get global vector x from DM  */
 76:   DMCreateGlobalVector(user.da, &x);

 78:   ini_bou(x, &user);
 79:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer);
 80:   VecView(x, viewer);
 81:   PetscViewerDestroy(&viewer);

 83:   TSCreate(PETSC_COMM_WORLD, &ts);
 84:   TSSetDM(ts, user.da);
 85:   TSSetProblemType(ts, TS_NONLINEAR);
 86:   TSSetType(ts, TSARKIMEX);
 87:   TSSetIFunction(ts, NULL, IFunction, &user);
 88:   /*  TSSetIJacobian(ts,NULL,NULL,IJacobian,&user);  */
 89:   TSSetApplicationContext(ts, &user);
 90:   TSSetTimeStep(ts, .005);
 91:   TSSetFromOptions(ts);
 92:   TSSetPostStep(ts, PostStep);
 93:   TSSolve(ts, x);

 95:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer);
 96:   VecView(x, viewer);
 97:   PetscViewerDestroy(&viewer);

 99:   VecDestroy(&x);
100:   DMDestroy(&user.da);
101:   TSDestroy(&ts);
102:   PetscFinalize();
103:   return 0;
104: }

106: PetscErrorCode PostStep(TS ts)
107: {
108:   Vec          X, gc;
109:   AppCtx      *user;
110:   PetscScalar  sum = 0, asum;
111:   PetscReal    t, **p;
112:   DMDACoor2d **coors;
113:   DM           cda;
114:   PetscInt     i, j, xs, ys, xm, ym;

116:   TSGetApplicationContext(ts, &user);
117:   TSGetTime(ts, &t);
118:   TSGetSolution(ts, &X);

120:   DMGetCoordinateDM(user->da, &cda);
121:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);
122:   DMGetCoordinates(user->da, &gc);
123:   DMDAVecGetArrayRead(cda, gc, &coors);
124:   DMDAVecGetArrayRead(user->da, X, &p);
125:   for (i = xs; i < xs + xm; i++) {
126:     for (j = ys; j < ys + ym; j++) {
127:       if (coors[j][i].y < 5) sum += p[j][i];
128:     }
129:   }
130:   DMDAVecRestoreArrayRead(cda, gc, &coors);
131:   DMDAVecRestoreArrayRead(user->da, X, &p);
132:   MPI_Allreduce(&sum, &asum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)ts));
133:   PetscPrintf(PETSC_COMM_WORLD, "sum(p)*dw*dtheta at t = %f = %f\n", (double)t, (double)(asum));
134:   if (sum < 1.0e-2) {
135:     TSSetConvergedReason(ts, TS_CONVERGED_USER);
136:     PetscPrintf(PETSC_COMM_WORLD, "Exiting TS as the integral of PDF is almost zero\n");
137:   }
138:   return 0;
139: }

141: PetscErrorCode ini_bou(Vec X, AppCtx *user)
142: {
143:   DM            cda;
144:   DMDACoor2d  **coors;
145:   PetscScalar **p;
146:   Vec           gc;
147:   PetscInt      i, j;
148:   PetscInt      xs, ys, xm, ym, M, N;
149:   PetscScalar   xi, yi;
150:   PetscScalar   sigmax = user->sigmax, sigmay = user->sigmay;
151:   PetscScalar   rho = user->rho;
152:   PetscScalar   muy = user->muy, mux;
153:   PetscMPIInt   rank;
154:   PetscScalar   sum;

157:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
158:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
159:   user->dx = (user->xmax - user->xmin) / (M - 1);
160:   user->dy = (user->ymax - user->ymin) / (N - 1);
161:   DMGetCoordinateDM(user->da, &cda);
162:   DMGetCoordinates(user->da, &gc);
163:   DMDAVecGetArray(cda, gc, &coors);
164:   DMDAVecGetArray(user->da, X, &p);
165:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

167:   /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable
168:      muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point
169:      in the y-direction. We only modify mux here
170:   */
171:   mux = user->mux = coors[0][M / 2 + 10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */
172:   if (user->nonoiseinitial) {
173:     for (i = xs; i < xs + xm; i++) {
174:       for (j = ys; j < ys + ym; j++) {
175:         xi = coors[j][i].x;
176:         yi = coors[j][i].y;
177:         if ((xi == mux) && (yi == muy)) p[j][i] = 1.0;
178:       }
179:     }
180:   } else {
181:     /* Change PM_min accordingly */
182:     user->PM_min = user->Pmax * PetscSinScalar(mux);
183:     for (i = xs; i < xs + xm; i++) {
184:       for (j = ys; j < ys + ym; j++) {
185:         xi = coors[j][i].x;
186:         yi = coors[j][i].y;
187:         p[j][i] = (0.5 / (PETSC_PI * sigmax * sigmay * PetscSqrtScalar(1.0 - rho * rho))) * PetscExpScalar(-0.5 / (1 - rho * rho) * (PetscPowScalar((xi - mux) / sigmax, 2) + PetscPowScalar((yi - muy) / sigmay, 2) - 2 * rho * (xi - mux) * (yi - muy) / (sigmax * sigmay)));
188:       }
189:     }
190:   }
191:   DMDAVecRestoreArray(cda, gc, &coors);
192:   DMDAVecRestoreArray(user->da, X, &p);
193:   VecSum(X, &sum);
194:   VecScale(X, 1.0 / sum);
195:   return 0;
196: }

198: /* First advection term */
199: PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
200: {
201:   PetscScalar f, fpos, fneg;
202:   f    = (y - user->ws);
203:   fpos = PetscMax(f, 0);
204:   fneg = PetscMin(f, 0);
205:   if (user->st_width == 1) {
206:     *p1 = fpos * (p[j][i] - p[j][i - 1]) / user->dx + fneg * (p[j][i + 1] - p[j][i]) / user->dx;
207:   } else if (user->st_width == 2) {
208:     *p1 = fpos * (3 * p[j][i] - 4 * p[j][i - 1] + p[j][i - 2]) / (2 * user->dx) + fneg * (-p[j][i + 2] + 4 * p[j][i + 1] - 3 * p[j][i]) / (2 * user->dx);
209:   } else if (user->st_width == 3) {
210:     *p1 = fpos * (2 * p[j][i + 1] + 3 * p[j][i] - 6 * p[j][i - 1] + p[j][i - 2]) / (6 * user->dx) + fneg * (-p[j][i + 2] + 6 * p[j][i + 1] - 3 * p[j][i] - 2 * p[j][i - 1]) / (6 * user->dx);
211:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
212:   return 0;
213: }

215: /* Second advection term */
216: PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
217: {
218:   PetscScalar f, fpos, fneg;
219:   f    = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x));
220:   fpos = PetscMax(f, 0);
221:   fneg = PetscMin(f, 0);
222:   if (user->st_width == 1) {
223:     *p2 = fpos * (p[j][i] - p[j - 1][i]) / user->dy + fneg * (p[j + 1][i] - p[j][i]) / user->dy;
224:   } else if (user->st_width == 2) {
225:     *p2 = fpos * (3 * p[j][i] - 4 * p[j - 1][i] + p[j - 2][i]) / (2 * user->dy) + fneg * (-p[j + 2][i] + 4 * p[j + 1][i] - 3 * p[j][i]) / (2 * user->dy);
226:   } else if (user->st_width == 3) {
227:     *p2 = fpos * (2 * p[j + 1][i] + 3 * p[j][i] - 6 * p[j - 1][i] + p[j - 2][i]) / (6 * user->dy) + fneg * (-p[j + 2][i] + 6 * p[j + 1][i] - 3 * p[j][i] - 2 * p[j - 1][i]) / (6 * user->dy);
228:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
229:   return 0;
230: }

232: /* Diffusion term */
233: PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
234: {
236:   if (user->st_width == 1) {
237:     *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
238:   } else if (user->st_width == 2) {
239:     *p_diff = user->disper_coe * ((-p[j - 2][i] + 16 * p[j - 1][i] - 30 * p[j][i] + 16 * p[j + 1][i] - p[j + 2][i]) / (12.0 * user->dy * user->dy));
240:   } else if (user->st_width == 3) {
241:     *p_diff = user->disper_coe * ((2 * p[j - 3][i] - 27 * p[j - 2][i] + 270 * p[j - 1][i] - 490 * p[j][i] + 270 * p[j + 1][i] - 27 * p[j + 2][i] + 2 * p[j + 3][i]) / (180.0 * user->dy * user->dy));
242:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
243:   return 0;
244: }

246: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
247: {
248:   AppCtx       *user = (AppCtx *)ctx;
249:   DM            cda;
250:   DMDACoor2d  **coors;
251:   PetscScalar **p, **f, **pdot;
252:   PetscInt      i, j;
253:   PetscInt      xs, ys, xm, ym, M, N;
254:   Vec           localX, gc, localXdot;
255:   PetscScalar   p_adv1 = 0.0, p_adv2 = 0.0, p_diff; /* initialize to prevent incorrect compiler warnings */

258:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
259:   DMGetCoordinateDM(user->da, &cda);
260:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

262:   DMGetLocalVector(user->da, &localX);
263:   DMGetLocalVector(user->da, &localXdot);

265:   DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX);
266:   DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX);
267:   DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot);
268:   DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot);

270:   DMGetCoordinatesLocal(user->da, &gc);

272:   DMDAVecGetArrayRead(cda, gc, &coors);
273:   DMDAVecGetArrayRead(user->da, localX, &p);
274:   DMDAVecGetArrayRead(user->da, localXdot, &pdot);
275:   DMDAVecGetArray(user->da, F, &f);

277:   user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda));
278:   for (i = xs; i < xs + xm; i++) {
279:     for (j = ys; j < ys + ym; j++) {
280:       adv1(p, coors[j][i].y, i, j, M, &p_adv1, user);
281:       adv2(p, coors[j][i].x, i, j, N, &p_adv2, user);
282:       diffuse(p, i, j, t, &p_diff, user);
283:       f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
284:     }
285:   }
286:   DMDAVecRestoreArrayRead(user->da, localX, &p);
287:   DMDAVecRestoreArrayRead(user->da, localX, &pdot);
288:   DMRestoreLocalVector(user->da, &localX);
289:   DMRestoreLocalVector(user->da, &localXdot);
290:   DMDAVecRestoreArray(user->da, F, &f);
291:   DMDAVecRestoreArrayRead(cda, gc, &coors);

293:   return 0;
294: }

296: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx)
297: {
298:   AppCtx      *user = (AppCtx *)ctx;
299:   DM           cda;
300:   DMDACoor2d **coors;
301:   PetscInt     i, j;
302:   PetscInt     xs, ys, xm, ym, M, N;
303:   Vec          gc;
304:   PetscScalar  val[5], xi, yi;
305:   MatStencil   row, col[5];
306:   PetscScalar  c1, c3, c5, c1pos, c1neg, c3pos, c3neg;

309:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
310:   DMGetCoordinateDM(user->da, &cda);
311:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

313:   DMGetCoordinatesLocal(user->da, &gc);
314:   DMDAVecGetArrayRead(cda, gc, &coors);
315:   for (i = xs; i < xs + xm; i++) {
316:     for (j = ys; j < ys + ym; j++) {
317:       PetscInt nc = 0;
318:       xi          = coors[j][i].x;
319:       yi          = coors[j][i].y;
320:       row.i       = i;
321:       row.j       = j;
322:       c1          = (yi - user->ws) / user->dx;
323:       c1pos       = PetscMax(c1, 0);
324:       c1neg       = PetscMin(c1, 0);
325:       c3          = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / user->dy;
326:       c3pos       = PetscMax(c3, 0);
327:       c3neg       = PetscMin(c3, 0);
328:       c5          = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
329:       col[nc].i   = i - 1;
330:       col[nc].j   = j;
331:       val[nc++]   = c1pos;
332:       col[nc].i   = i + 1;
333:       col[nc].j   = j;
334:       val[nc++]   = -c1neg;
335:       col[nc].i   = i;
336:       col[nc].j   = j - 1;
337:       val[nc++]   = c3pos + c5;
338:       col[nc].i   = i;
339:       col[nc].j   = j + 1;
340:       val[nc++]   = -c3neg + c5;
341:       col[nc].i   = i;
342:       col[nc].j   = j;
343:       val[nc++]   = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a;
344:       MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES);
345:     }
346:   }
347:   DMDAVecRestoreArrayRead(cda, gc, &coors);

349:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
350:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
351:   if (J != Jpre) {
352:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
353:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
354:   }
355:   return 0;
356: }

358: PetscErrorCode Parameter_settings(AppCtx *user)
359: {
360:   PetscBool flg;


364:   /* Set default parameters */
365:   user->ws             = 1.0;
366:   user->H              = 5.0;
367:   user->Pmax           = 2.1;
368:   user->PM_min         = 1.0;
369:   user->lambda         = 0.1;
370:   user->q              = 1.0;
371:   user->mux            = PetscAsinScalar(user->PM_min / user->Pmax);
372:   user->sigmax         = 0.1;
373:   user->sigmay         = 0.1;
374:   user->rho            = 0.0;
375:   user->xmin           = -PETSC_PI;
376:   user->xmax           = PETSC_PI;
377:   user->bx             = DM_BOUNDARY_PERIODIC;
378:   user->by             = DM_BOUNDARY_MIRROR;
379:   user->nonoiseinitial = PETSC_FALSE;

381:   /*
382:      ymin of -3 seems to let the unstable solution move up and leave a zero in its wake
383:      with an ymin of -1 the wake is never exactly zero
384:   */
385:   user->ymin     = -3.0;
386:   user->ymax     = 10.0;
387:   user->st_width = 1;

389:   PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg);
390:   PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg);
391:   PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg);
392:   PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg);
393:   PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg);
394:   PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg);
395:   PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg);
396:   PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg);
397:   PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg);
398:   if (flg == 0) user->muy = user->ws;
399:   PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg);
400:   PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg);
401:   PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg);
402:   PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg);
403:   PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg);
404:   PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg);
405:   PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg);
406:   PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg);
407:   PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg);
408:   PetscOptionsGetBool(NULL, NULL, "-nonoiseinitial", &user->nonoiseinitial, &flg);

410:   return 0;
411: }

413: /*TEST

415:    build:
416:       requires: !complex !single

418:    test:
419:       args: -ts_max_steps 2
420:       localrunfiles: petscopt_ex7

422:    test:
423:       suffix: 2
424:       args: -ts_max_steps 2 -snes_mf_operator
425:       output_file: output/ex7_1.out
426:       localrunfiles: petscopt_ex7
427:       timeoutfactor: 2

429:    test:
430:       suffix: 3
431:       args: -ts_max_steps 2 -snes_mf -pc_type none
432:       output_file: output/ex7_1.out
433:       localrunfiles: petscopt_ex7
434:       timeoutfactor: 2

436: TEST*/