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;
5
6 Boundary conditions:
7 Zero dirichlet in y using ghosted values
8 Periodic in x
9
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))
13
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.
16
17 */
18
19 #include <petscdm.h>
20 #include <petscdmda.h>
21 #include <petscts.h>
22
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;
52
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);
58
main(int argc,char ** argv)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;
65
66 PetscFunctionBeginUser;
67 PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex8", help));
68
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"));
79
80 /* Get global vector x from DM */
81 PetscCall(DMCreateGlobalVector(user.da, &x));
82
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));
87
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));
99
100 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer));
101 PetscCall(VecView(x, viewer));
102 PetscCall(PetscViewerDestroy(&viewer));
103
104 PetscCall(VecDestroy(&x));
105 PetscCall(DMDestroy(&user.da));
106 PetscCall(TSDestroy(&ts));
107 PetscCall(PetscFinalize());
108 return 0;
109 }
110
PostStep(TS ts)111 PetscErrorCode PostStep(TS ts)
112 {
113 Vec X;
114 AppCtx *user;
115 PetscReal t;
116 PetscScalar asum;
117
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 }*/
129
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;
132
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 }
137
ini_bou(Vec X,AppCtx * user)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;
146
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));
156
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;
164
165 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
166 PetscCall(DMDAVecRestoreArray(user->da, X, &p));
167 PetscFunctionReturn(PETSC_SUCCESS);
168 }
169
170 /* First advection term */
adv1(PetscScalar ** p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar * p1,AppCtx * user)171 PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
172 {
173 PetscScalar f = (y - user->ws), fpos = PetscMax(f, 0), fneg = PetscMin(f, 0);
174
175 PetscFunctionBegin;
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 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
187 /* Second advection term */
adv2(PetscScalar ** p,PetscScalar x,PetscScalar y,PetscInt i,PetscInt j,PetscInt N,PetscScalar * p2,AppCtx * user)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
192 PetscFunctionBegin;
193 f = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x) - user->D * (y - user->ws));
194 fpos = PetscMax(f, 0);
195 fneg = PetscMin(f, 0);
196 if (user->st_width == 1) {
197 *p2 = fpos * (p[j][i] - p[j - 1][i]) / user->dy + fneg * (p[j + 1][i] - p[j][i]) / user->dy;
198 } else if (user->st_width == 2) {
199 *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);
200 } else if (user->st_width == 3) {
201 *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);
202 }
203
204 /* *p2 = f*(p[j+1][i] - p[j-1][i])/user->dy;*/
205 PetscFunctionReturn(PETSC_SUCCESS);
206 }
207
208 /* Diffusion term */
diffuse(PetscScalar ** p,PetscInt i,PetscInt j,PetscReal t,PetscScalar * p_diff,AppCtx * user)209 PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
210 {
211 PetscFunctionBeginUser;
212 if (user->st_width == 1) {
213 *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
214 } else if (user->st_width == 2) {
215 *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));
216 } else if (user->st_width == 3) {
217 *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));
218 }
219 PetscFunctionReturn(PETSC_SUCCESS);
220 }
221
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,PetscCtx ctx)222 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, PetscCtx ctx)
223 {
224 AppCtx *user = (AppCtx *)ctx;
225 DM cda;
226 DMDACoor2d **coors;
227 PetscScalar **p, **f, **pdot;
228 PetscInt i, j;
229 PetscInt xs, ys, xm, ym, M, N;
230 Vec localX, gc, localXdot;
231 PetscScalar p_adv1 = 0.0, p_adv2 = 0.0, p_diff = 0;
232 PetscScalar diffuse1, gamma;
233
234 PetscFunctionBeginUser;
235 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
236 PetscCall(DMGetCoordinateDM(user->da, &cda));
237 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
238
239 PetscCall(DMGetLocalVector(user->da, &localX));
240 PetscCall(DMGetLocalVector(user->da, &localXdot));
241
242 PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX));
243 PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX));
244 PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot));
245 PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot));
246
247 PetscCall(DMGetCoordinatesLocal(user->da, &gc));
248
249 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
250 PetscCall(DMDAVecGetArrayRead(user->da, localX, &p));
251 PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot));
252 PetscCall(DMDAVecGetArray(user->da, F, &f));
253
254 gamma = user->D * user->ws / (2 * user->H);
255 diffuse1 = user->lambda * user->lambda * user->q / (user->lambda * gamma + 1) * (1.0 - PetscExpScalar(-t * (gamma + 1.0) / user->lambda));
256 user->disper_coe = user->ws * user->ws / (4 * user->H * user->H) * diffuse1;
257
258 for (i = xs; i < xs + xm; i++) {
259 for (j = ys; j < ys + ym; j++) {
260 PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user));
261 PetscCall(adv2(p, coors[j][i].x, coors[j][i].y, i, j, N, &p_adv2, user));
262 PetscCall(diffuse(p, i, j, t, &p_diff, user));
263 f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
264 }
265 }
266 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p));
267 PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot));
268 PetscCall(DMRestoreLocalVector(user->da, &localX));
269 PetscCall(DMRestoreLocalVector(user->da, &localXdot));
270 PetscCall(DMDAVecRestoreArray(user->da, F, &f));
271 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
272 PetscFunctionReturn(PETSC_SUCCESS);
273 }
274
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,PetscCtx ctx)275 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, PetscCtx ctx)
276 {
277 AppCtx *user = (AppCtx *)ctx;
278 DM cda;
279 DMDACoor2d **coors;
280 PetscInt i, j;
281 PetscInt xs, ys, xm, ym, M, N;
282 Vec gc;
283 PetscScalar val[5], xi, yi;
284 MatStencil row, col[5];
285 PetscScalar c1, c3, c5, c1pos, c1neg, c3pos, c3neg;
286
287 PetscFunctionBeginUser;
288 PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
289 PetscCall(DMGetCoordinateDM(user->da, &cda));
290 PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
291
292 PetscCall(DMGetCoordinatesLocal(user->da, &gc));
293 PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
294 for (i = xs; i < xs + xm; i++) {
295 for (j = ys; j < ys + ym; j++) {
296 PetscInt nc = 0;
297 xi = coors[j][i].x;
298 yi = coors[j][i].y;
299 row.i = i;
300 row.j = j;
301 c1 = (yi - user->ws) / user->dx;
302 c1pos = PetscMax(c1, 0);
303 c1neg = PetscMin(c1, 0);
304 c3 = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi) - user->D * (yi - user->ws)) / user->dy;
305 c3pos = PetscMax(c3, 0);
306 c3neg = PetscMin(c3, 0);
307 c5 = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
308 col[nc].i = i - 1;
309 col[nc].j = j;
310 val[nc++] = c1pos;
311 col[nc].i = i + 1;
312 col[nc].j = j;
313 val[nc++] = -c1neg;
314 col[nc].i = i;
315 col[nc].j = j - 1;
316 val[nc++] = c3pos + c5;
317 col[nc].i = i;
318 col[nc].j = j + 1;
319 val[nc++] = -c3neg + c5;
320 col[nc].i = i;
321 col[nc].j = j;
322 val[nc++] = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a;
323 PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
324 }
325 }
326 PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
327
328 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
329 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
330 if (J != Jpre) {
331 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
332 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
333 }
334 PetscFunctionReturn(PETSC_SUCCESS);
335 }
336
Parameter_settings(AppCtx * user)337 PetscErrorCode Parameter_settings(AppCtx *user)
338 {
339 PetscBool flg;
340
341 PetscFunctionBeginUser;
342 /* Set default parameters */
343 user->ws = 1.0;
344 user->H = 5.0;
345 user->D = 0.0;
346 user->Pmax = user->Pmax_s = 2.1;
347 user->PM_min = 1.0;
348 user->lambda = 0.1;
349 user->q = 1.0;
350 user->mux = PetscAsinScalar(user->PM_min / user->Pmax);
351 user->sigmax = 0.1;
352 user->sigmay = 0.1;
353 user->rho = 0.0;
354 user->xmin = -PETSC_PI;
355 user->xmax = PETSC_PI;
356 user->bx = DM_BOUNDARY_PERIODIC;
357 user->by = DM_BOUNDARY_GHOSTED;
358 user->tf = user->tcl = -1;
359 user->ymin = -2.0;
360 user->ymax = 2.0;
361 user->st_width = 1;
362
363 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg));
364 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg));
365 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-D", &user->D, &flg));
366 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg));
367 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg));
368 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg));
369 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg));
370 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg));
371 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg));
372 if (flg == 0) user->muy = user->ws;
373 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg));
374 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg));
375 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg));
376 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg));
377 PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg));
378 PetscCall(PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg));
379 PetscCall(PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg));
380 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tf", &user->tf, &flg));
381 PetscCall(PetscOptionsGetReal(NULL, NULL, "-tcl", &user->tcl, &flg));
382 PetscFunctionReturn(PETSC_SUCCESS);
383 }
384
385 /*TEST
386
387 build:
388 requires: !complex x
389
390 test:
391 args: -ts_max_steps 1
392 localrunfiles: petscopt_ex8
393 timeoutfactor: 3
394
395 TEST*/
396