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;
66: PetscFunctionBeginUser;
67: PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex8", help));
69: /* Get physics and time parameters */
70: PetscCall(Parameter_settings(&user));
71: /* Create a 2D DA with dof = 1 */
72: PetscCall(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: PetscCall(DMSetFromOptions(user.da));
74: PetscCall(DMSetUp(user.da));
75: /* Set x and y coordinates */
76: PetscCall(DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0, 0));
77: PetscCall(DMDASetCoordinateName(user.da, 0, "X - the angle"));
78: PetscCall(DMDASetCoordinateName(user.da, 1, "Y - the speed"));
80: /* Get global vector x from DM */
81: PetscCall(DMCreateGlobalVector(user.da, &x));
83: PetscCall(ini_bou(x, &user));
84: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer));
85: PetscCall(VecView(x, viewer));
86: PetscCall(PetscViewerDestroy(&viewer));
88: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
89: PetscCall(TSSetDM(ts, user.da));
90: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
91: PetscCall(TSSetType(ts, TSARKIMEX));
92: PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
93: /* PetscCall(TSSetIJacobian(ts,NULL,NULL,IJacobian,&user)); */
94: PetscCall(TSSetApplicationContext(ts, &user));
95: PetscCall(TSSetTimeStep(ts, .005));
96: PetscCall(TSSetFromOptions(ts));
97: PetscCall(TSSetPostStep(ts, PostStep));
98: PetscCall(TSSolve(ts, x));
100: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer));
101: PetscCall(VecView(x, viewer));
102: PetscCall(PetscViewerDestroy(&viewer));
104: PetscCall(VecDestroy(&x));
105: PetscCall(DMDestroy(&user.da));
106: PetscCall(TSDestroy(&ts));
107: PetscCall(PetscFinalize());
108: return 0;
109: }
111: PetscErrorCode PostStep(TS ts)
112: {
113: Vec X;
114: AppCtx *user;
115: PetscReal t;
116: PetscScalar asum;
118: PetscFunctionBegin;
119: PetscCall(TSGetApplicationContext(ts, &user));
120: PetscCall(TSGetTime(ts, &t));
121: PetscCall(TSGetSolution(ts, &X));
122: /*
123: if (t >= .2) {
124: PetscCall(TSGetSolution(ts,&X));
125: PetscCall(VecView(X,PETSC_VIEWER_BINARY_WORLD));
126: exit(0);
127: results in initial conditions after fault in binaryoutput
128: }*/
130: 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 */
131: else user->Pmax = user->Pmax_s;
133: PetscCall(VecSum(X, &asum));
134: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sum(p) at t = %f = %f\n", (double)t, (double)(asum)));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: PetscErrorCode ini_bou(Vec X, AppCtx *user)
139: {
140: DM cda;
141: DMDACoor2d **coors;
142: PetscScalar **p;
143: Vec gc;
144: PetscInt M, N, Ir, J;
145: PetscMPIInt rank;
147: PetscFunctionBeginUser;
148: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
149: PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
150: user->dx = (user->xmax - user->xmin) / (M - 1);
151: user->dy = (user->ymax - user->ymin) / (N - 1);
152: PetscCall(DMGetCoordinateDM(user->da, &cda));
153: PetscCall(DMGetCoordinates(user->da, &gc));
154: PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
155: PetscCall(DMDAVecGetArray(user->da, X, &p));
157: /* Point mass at (mux,muy) */
158: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original user->mux = %f, user->muy = %f\n", (double)PetscRealPart(user->mux), (double)PetscRealPart(user->muy)));
159: PetscCall(DMDAGetLogicalCoordinate(user->da, user->mux, user->muy, 0.0, &Ir, &J, NULL, &user->mux, &user->muy, NULL));
160: user->PM_min = user->Pmax * PetscSinScalar(user->mux);
161: PetscCall(
162: 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)));
163: if (Ir > -1 && J > -1) p[J][Ir] = 1.0;
165: PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
166: PetscCall(DMDAVecRestoreArray(user->da, X, &p));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /* First advection term */
171: PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
172: {
173: PetscScalar f, fpos, fneg;
174: PetscFunctionBegin;
175: f = (y - user->ws);
176: fpos = PetscMax(f, 0);
177: fneg = PetscMin(f, 0);
178: if (user->st_width == 1) {
179: *p1 = fpos * (p[j][i] - p[j][i - 1]) / user->dx + fneg * (p[j][i + 1] - p[j][i]) / user->dx;
180: } else if (user->st_width == 2) {
181: *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);
182: } else if (user->st_width == 3) {
183: *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);
184: }
185: /* *p1 = f*(p[j][i+1] - p[j][i-1])/user->dx;*/
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /* Second advection term */
190: PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscScalar y, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
191: {
192: PetscScalar f, fpos, fneg;
193: PetscFunctionBegin;
194: f = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x) - user->D * (y - user->ws));
195: fpos = PetscMax(f, 0);
196: fneg = PetscMin(f, 0);
197: if (user->st_width == 1) {
198: *p2 = fpos * (p[j][i] - p[j - 1][i]) / user->dy + fneg * (p[j + 1][i] - p[j][i]) / user->dy;
199: } else if (user->st_width == 2) {
200: *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);
201: } else if (user->st_width == 3) {
202: *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);
203: }
205: /* *p2 = f*(p[j+1][i] - p[j-1][i])/user->dy;*/
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /* Diffusion term */
210: PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
211: {
212: PetscFunctionBeginUser;
213: if (user->st_width == 1) {
214: *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
215: } else if (user->st_width == 2) {
216: *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));
217: } else if (user->st_width == 3) {
218: *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));
219: }
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
224: {
225: AppCtx *user = (AppCtx *)ctx;
226: DM cda;
227: DMDACoor2d **coors;
228: PetscScalar **p, **f, **pdot;
229: PetscInt i, j;
230: PetscInt xs, ys, xm, ym, M, N;
231: Vec localX, gc, localXdot;
232: PetscScalar p_adv1 = 0.0, p_adv2 = 0.0, p_diff = 0;
233: PetscScalar diffuse1, gamma;
235: PetscFunctionBeginUser;
236: PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
237: PetscCall(DMGetCoordinateDM(user->da, &cda));
238: PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
240: PetscCall(DMGetLocalVector(user->da, &localX));
241: PetscCall(DMGetLocalVector(user->da, &localXdot));
243: PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX));
244: PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX));
245: PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot));
246: PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot));
248: PetscCall(DMGetCoordinatesLocal(user->da, &gc));
250: PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
251: PetscCall(DMDAVecGetArrayRead(user->da, localX, &p));
252: PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot));
253: PetscCall(DMDAVecGetArray(user->da, F, &f));
255: gamma = user->D * user->ws / (2 * user->H);
256: diffuse1 = user->lambda * user->lambda * user->q / (user->lambda * gamma + 1) * (1.0 - PetscExpScalar(-t * (gamma + 1.0) / user->lambda));
257: user->disper_coe = user->ws * user->ws / (4 * user->H * user->H) * diffuse1;
259: for (i = xs; i < xs + xm; i++) {
260: for (j = ys; j < ys + ym; j++) {
261: PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user));
262: PetscCall(adv2(p, coors[j][i].x, coors[j][i].y, i, j, N, &p_adv2, user));
263: PetscCall(diffuse(p, i, j, t, &p_diff, user));
264: f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
265: }
266: }
267: PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p));
268: PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot));
269: PetscCall(DMRestoreLocalVector(user->da, &localX));
270: PetscCall(DMRestoreLocalVector(user->da, &localXdot));
271: PetscCall(DMDAVecRestoreArray(user->da, F, &f));
272: PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ctx)
278: {
279: AppCtx *user = (AppCtx *)ctx;
280: DM cda;
281: DMDACoor2d **coors;
282: PetscInt i, j;
283: PetscInt xs, ys, xm, ym, M, N;
284: Vec gc;
285: PetscScalar val[5], xi, yi;
286: MatStencil row, col[5];
287: PetscScalar c1, c3, c5, c1pos, c1neg, c3pos, c3neg;
289: PetscFunctionBeginUser;
290: PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
291: PetscCall(DMGetCoordinateDM(user->da, &cda));
292: PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
294: PetscCall(DMGetCoordinatesLocal(user->da, &gc));
295: PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
296: for (i = xs; i < xs + xm; i++) {
297: for (j = ys; j < ys + ym; j++) {
298: PetscInt nc = 0;
299: xi = coors[j][i].x;
300: yi = coors[j][i].y;
301: row.i = i;
302: row.j = j;
303: c1 = (yi - user->ws) / user->dx;
304: c1pos = PetscMax(c1, 0);
305: c1neg = PetscMin(c1, 0);
306: c3 = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi) - user->D * (yi - user->ws)) / user->dy;
307: c3pos = PetscMax(c3, 0);
308: c3neg = PetscMin(c3, 0);
309: c5 = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
310: col[nc].i = i - 1;
311: col[nc].j = j;
312: val[nc++] = c1pos;
313: col[nc].i = i + 1;
314: col[nc].j = j;
315: val[nc++] = -c1neg;
316: col[nc].i = i;
317: col[nc].j = j - 1;
318: val[nc++] = c3pos + c5;
319: col[nc].i = i;
320: col[nc].j = j + 1;
321: val[nc++] = -c3neg + c5;
322: col[nc].i = i;
323: col[nc].j = j;
324: val[nc++] = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a;
325: PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
326: }
327: }
328: PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
330: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
331: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
332: if (J != Jpre) {
333: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
334: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
335: }
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PetscErrorCode Parameter_settings(AppCtx *user)
340: {
341: PetscBool flg;
343: PetscFunctionBeginUser;
344: /* Set default parameters */
345: user->ws = 1.0;
346: user->H = 5.0;
347: user->D = 0.0;
348: user->Pmax = user->Pmax_s = 2.1;
349: user->PM_min = 1.0;
350: user->lambda = 0.1;
351: user->q = 1.0;
352: user->mux = PetscAsinScalar(user->PM_min / user->Pmax);
353: user->sigmax = 0.1;
354: user->sigmay = 0.1;
355: user->rho = 0.0;
356: user->xmin = -PETSC_PI;
357: user->xmax = PETSC_PI;
358: user->bx = DM_BOUNDARY_PERIODIC;
359: user->by = DM_BOUNDARY_GHOSTED;
360: user->tf = user->tcl = -1;
361: user->ymin = -2.0;
362: user->ymax = 2.0;
363: user->st_width = 1;
365: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg));
366: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg));
367: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-D", &user->D, &flg));
368: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg));
369: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg));
370: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg));
371: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg));
372: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg));
373: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg));
374: if (flg == 0) user->muy = user->ws;
375: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg));
376: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg));
377: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg));
378: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg));
379: PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg));
380: PetscCall(PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg));
381: PetscCall(PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg));
382: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tf", &user->tf, &flg));
383: PetscCall(PetscOptionsGetReal(NULL, NULL, "-tcl", &user->tcl, &flg));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*TEST
389: build:
390: requires: !complex x
392: test:
393: args: -ts_max_steps 1
394: localrunfiles: petscopt_ex8
395: timeoutfactor: 3
397: TEST*/