Actual source code: ex8.c

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

  6:    Boundary conditions:
  7:    Zero dirichlet in y using ghosted values
  8:    Periodic in x

 10:    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.
 11:    x_t = (y - ws)
 12:    y_t = (ws/2H)*(Pm - Pmax*sin(x) - D*(w - ws))

 14:    In this example, we can see the effect of a fault, that zeroes the electrical power output
 15:    Pmax*sin(x), on the PDF. The fault on/off times can be controlled by options -tf and -tcl respectively.

 17: */

 19: #include <petscdm.h>
 20: #include <petscdmda.h>
 21: #include <petscts.h>

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

 53: PetscErrorCode Parameter_settings(AppCtx *);
 54: PetscErrorCode ini_bou(Vec, AppCtx *);
 55: PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
 56: PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
 57: PetscErrorCode PostStep(TS);

 59: int main(int argc, char **argv)
 60: {
 61:   Vec         x;    /* Solution vector */
 62:   TS          ts;   /* Time-stepping context */
 63:   AppCtx      user; /* Application context */
 64:   PetscViewer viewer;

 67:   PetscInitialize(&argc, &argv, "petscopt_ex8", help);

 69:   /* Get physics and time parameters */
 70:   Parameter_settings(&user);
 71:   /* Create a 2D DA with dof = 1 */
 72:   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);
 73:   DMSetFromOptions(user.da);
 74:   DMSetUp(user.da);
 75:   /* Set x and y coordinates */
 76:   DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0, 0);
 77:   DMDASetCoordinateName(user.da, 0, "X - the angle");
 78:   DMDASetCoordinateName(user.da, 1, "Y - the speed");

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

 83:   ini_bou(x, &user);
 84:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer);
 85:   VecView(x, viewer);
 86:   PetscViewerDestroy(&viewer);

 88:   TSCreate(PETSC_COMM_WORLD, &ts);
 89:   TSSetDM(ts, user.da);
 90:   TSSetProblemType(ts, TS_NONLINEAR);
 91:   TSSetType(ts, TSARKIMEX);
 92:   TSSetIFunction(ts, NULL, IFunction, &user);
 93:   /*  TSSetIJacobian(ts,NULL,NULL,IJacobian,&user); */
 94:   TSSetApplicationContext(ts, &user);
 95:   TSSetTimeStep(ts, .005);
 96:   TSSetFromOptions(ts);
 97:   TSSetPostStep(ts, PostStep);
 98:   TSSolve(ts, x);

100:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer);
101:   VecView(x, viewer);
102:   PetscViewerDestroy(&viewer);

104:   VecDestroy(&x);
105:   DMDestroy(&user.da);
106:   TSDestroy(&ts);
107:   PetscFinalize();
108:   return 0;
109: }

111: PetscErrorCode PostStep(TS ts)
112: {
113:   Vec         X;
114:   AppCtx     *user;
115:   PetscReal   t;
116:   PetscScalar asum;

118:   TSGetApplicationContext(ts, &user);
119:   TSGetTime(ts, &t);
120:   TSGetSolution(ts, &X);
121:   /*
122:   if (t >= .2) {
123:     TSGetSolution(ts,&X);
124:     VecView(X,PETSC_VIEWER_BINARY_WORLD);
125:     exit(0);
126:      results in initial conditions after fault in binaryoutput
127:   }*/

129:   if ((t > user->tf) && (t < user->tcl)) user->Pmax = 0.0; /* A short-circuit that drives the electrical power output (Pmax*sin(delta)) to zero */
130:   else user->Pmax = user->Pmax_s;

132:   VecSum(X, &asum);
133:   PetscPrintf(PETSC_COMM_WORLD, "sum(p) at t = %f = %f\n", (double)t, (double)(asum));
134:   return 0;
135: }

137: PetscErrorCode ini_bou(Vec X, AppCtx *user)
138: {
139:   DM            cda;
140:   DMDACoor2d  **coors;
141:   PetscScalar **p;
142:   Vec           gc;
143:   PetscInt      M, N, Ir, J;
144:   PetscMPIInt   rank;

147:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
148:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
149:   user->dx = (user->xmax - user->xmin) / (M - 1);
150:   user->dy = (user->ymax - user->ymin) / (N - 1);
151:   DMGetCoordinateDM(user->da, &cda);
152:   DMGetCoordinates(user->da, &gc);
153:   DMDAVecGetArrayRead(cda, gc, &coors);
154:   DMDAVecGetArray(user->da, X, &p);

156:   /* Point mass at (mux,muy) */
157:   PetscPrintf(PETSC_COMM_WORLD, "Original user->mux = %f, user->muy = %f\n", (double)PetscRealPart(user->mux), (double)PetscRealPart(user->muy));
158:   DMDAGetLogicalCoordinate(user->da, user->mux, user->muy, 0.0, &Ir, &J, NULL, &user->mux, &user->muy, NULL);
159:   user->PM_min = user->Pmax * PetscSinScalar(user->mux);
160:   PetscCall(
161:     PetscPrintf(PETSC_COMM_WORLD, "Corrected user->mux = %f, user->muy = %f user->PM_min = %f,user->dx = %f\n", (double)PetscRealPart(user->mux), (double)PetscRealPart(user->muy), (double)PetscRealPart(user->PM_min), (double)PetscRealPart(user->dx)));
162:   if (Ir > -1 && J > -1) p[J][Ir] = 1.0;

164:   DMDAVecRestoreArrayRead(cda, gc, &coors);
165:   DMDAVecRestoreArray(user->da, X, &p);
166:   return 0;
167: }

169: /* First advection term */
170: PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
171: {
172:   PetscScalar f, fpos, fneg;
173:   f    = (y - user->ws);
174:   fpos = PetscMax(f, 0);
175:   fneg = PetscMin(f, 0);
176:   if (user->st_width == 1) {
177:     *p1 = fpos * (p[j][i] - p[j][i - 1]) / user->dx + fneg * (p[j][i + 1] - p[j][i]) / user->dx;
178:   } else if (user->st_width == 2) {
179:     *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);
180:   } else if (user->st_width == 3) {
181:     *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);
182:   }
183:   /* *p1 = f*(p[j][i+1] - p[j][i-1])/user->dx;*/
184:   return 0;
185: }

187: /* Second advection term */
188: PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscScalar y, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
189: {
190:   PetscScalar f, fpos, fneg;
191:   f    = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x) - user->D * (y - user->ws));
192:   fpos = PetscMax(f, 0);
193:   fneg = PetscMin(f, 0);
194:   if (user->st_width == 1) {
195:     *p2 = fpos * (p[j][i] - p[j - 1][i]) / user->dy + fneg * (p[j + 1][i] - p[j][i]) / user->dy;
196:   } else if (user->st_width == 2) {
197:     *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);
198:   } else if (user->st_width == 3) {
199:     *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);
200:   }

202:   /* *p2 = f*(p[j+1][i] - p[j-1][i])/user->dy;*/
203:   return 0;
204: }

206: /* Diffusion term */
207: PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
208: {
210:   if (user->st_width == 1) {
211:     *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
212:   } else if (user->st_width == 2) {
213:     *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));
214:   } else if (user->st_width == 3) {
215:     *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));
216:   }
217:   return 0;
218: }

220: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
221: {
222:   AppCtx       *user = (AppCtx *)ctx;
223:   DM            cda;
224:   DMDACoor2d  **coors;
225:   PetscScalar **p, **f, **pdot;
226:   PetscInt      i, j;
227:   PetscInt      xs, ys, xm, ym, M, N;
228:   Vec           localX, gc, localXdot;
229:   PetscScalar   p_adv1 = 0.0, p_adv2 = 0.0, p_diff = 0;
230:   PetscScalar   diffuse1, gamma;

233:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
234:   DMGetCoordinateDM(user->da, &cda);
235:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

237:   DMGetLocalVector(user->da, &localX);
238:   DMGetLocalVector(user->da, &localXdot);

240:   DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX);
241:   DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX);
242:   DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot);
243:   DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot);

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

247:   DMDAVecGetArrayRead(cda, gc, &coors);
248:   DMDAVecGetArrayRead(user->da, localX, &p);
249:   DMDAVecGetArrayRead(user->da, localXdot, &pdot);
250:   DMDAVecGetArray(user->da, F, &f);

252:   gamma            = user->D * user->ws / (2 * user->H);
253:   diffuse1         = user->lambda * user->lambda * user->q / (user->lambda * gamma + 1) * (1.0 - PetscExpScalar(-t * (gamma + 1.0) / user->lambda));
254:   user->disper_coe = user->ws * user->ws / (4 * user->H * user->H) * diffuse1;

256:   for (i = xs; i < xs + xm; i++) {
257:     for (j = ys; j < ys + ym; j++) {
258:       adv1(p, coors[j][i].y, i, j, M, &p_adv1, user);
259:       adv2(p, coors[j][i].x, coors[j][i].y, i, j, N, &p_adv2, user);
260:       diffuse(p, i, j, t, &p_diff, user);
261:       f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
262:     }
263:   }
264:   DMDAVecRestoreArrayRead(user->da, localX, &p);
265:   DMDAVecRestoreArrayRead(user->da, localX, &pdot);
266:   DMRestoreLocalVector(user->da, &localX);
267:   DMRestoreLocalVector(user->da, &localXdot);
268:   DMDAVecRestoreArray(user->da, F, &f);
269:   DMDAVecRestoreArrayRead(cda, gc, &coors);

271:   return 0;
272: }

274: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx)
275: {
276:   AppCtx      *user = (AppCtx *)ctx;
277:   DM           cda;
278:   DMDACoor2d **coors;
279:   PetscInt     i, j;
280:   PetscInt     xs, ys, xm, ym, M, N;
281:   Vec          gc;
282:   PetscScalar  val[5], xi, yi;
283:   MatStencil   row, col[5];
284:   PetscScalar  c1, c3, c5, c1pos, c1neg, c3pos, c3neg;

287:   DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
288:   DMGetCoordinateDM(user->da, &cda);
289:   DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0);

291:   DMGetCoordinatesLocal(user->da, &gc);
292:   DMDAVecGetArrayRead(cda, gc, &coors);
293:   for (i = xs; i < xs + xm; i++) {
294:     for (j = ys; j < ys + ym; j++) {
295:       PetscInt nc = 0;
296:       xi          = coors[j][i].x;
297:       yi          = coors[j][i].y;
298:       row.i       = i;
299:       row.j       = j;
300:       c1          = (yi - user->ws) / user->dx;
301:       c1pos       = PetscMax(c1, 0);
302:       c1neg       = PetscMin(c1, 0);
303:       c3          = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi) - user->D * (yi - user->ws)) / user->dy;
304:       c3pos       = PetscMax(c3, 0);
305:       c3neg       = PetscMin(c3, 0);
306:       c5          = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
307:       col[nc].i   = i - 1;
308:       col[nc].j   = j;
309:       val[nc++]   = c1pos;
310:       col[nc].i   = i + 1;
311:       col[nc].j   = j;
312:       val[nc++]   = -c1neg;
313:       col[nc].i   = i;
314:       col[nc].j   = j - 1;
315:       val[nc++]   = c3pos + c5;
316:       col[nc].i   = i;
317:       col[nc].j   = j + 1;
318:       val[nc++]   = -c3neg + c5;
319:       col[nc].i   = i;
320:       col[nc].j   = j;
321:       val[nc++]   = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a;
322:       MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES);
323:     }
324:   }
325:   DMDAVecRestoreArrayRead(cda, gc, &coors);

327:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
328:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
329:   if (J != Jpre) {
330:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
331:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
332:   }
333:   return 0;
334: }

336: PetscErrorCode Parameter_settings(AppCtx *user)
337: {
338:   PetscBool flg;

341:   /* Set default parameters */
342:   user->ws   = 1.0;
343:   user->H    = 5.0;
344:   user->D    = 0.0;
345:   user->Pmax = user->Pmax_s = 2.1;
346:   user->PM_min              = 1.0;
347:   user->lambda              = 0.1;
348:   user->q                   = 1.0;
349:   user->mux                 = PetscAsinScalar(user->PM_min / user->Pmax);
350:   user->sigmax              = 0.1;
351:   user->sigmay              = 0.1;
352:   user->rho                 = 0.0;
353:   user->xmin                = -PETSC_PI;
354:   user->xmax                = PETSC_PI;
355:   user->bx                  = DM_BOUNDARY_PERIODIC;
356:   user->by                  = DM_BOUNDARY_GHOSTED;
357:   user->tf = user->tcl = -1;
358:   user->ymin           = -2.0;
359:   user->ymax           = 2.0;
360:   user->st_width       = 1;

362:   PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg);
363:   PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg);
364:   PetscOptionsGetScalar(NULL, NULL, "-D", &user->D, &flg);
365:   PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg);
366:   PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg);
367:   PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg);
368:   PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg);
369:   PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg);
370:   PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg);
371:   if (flg == 0) user->muy = user->ws;
372:   PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg);
373:   PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg);
374:   PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg);
375:   PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg);
376:   PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg);
377:   PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg);
378:   PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg);
379:   PetscOptionsGetReal(NULL, NULL, "-tf", &user->tf, &flg);
380:   PetscOptionsGetReal(NULL, NULL, "-tcl", &user->tcl, &flg);
381:   return 0;
382: }

384: /*TEST

386:    build:
387:       requires: !complex x

389:    test:
390:       args: -ts_max_steps 1
391:       localrunfiles: petscopt_ex8
392:       timeoutfactor: 3

394: TEST*/