Actual source code: ex7.c
petsc-3.12.2 2019-11-22
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: {
57: Vec x; /* Solution vector */
58: TS ts; /* Time-stepping context */
59: AppCtx user; /* Application context */
60: PetscViewer viewer;
62: PetscInitialize(&argc,&argv,"petscopt_ex7", help);if (ierr) return ierr;
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 ierr;
104: }
106: PetscErrorCode PostStep(TS ts)
107: {
109: Vec X,gc;
110: AppCtx *user;
111: PetscScalar sum = 0,asum;
112: PetscReal t,**p;
113: DMDACoor2d **coors;
114: DM cda;
115: PetscInt i,j,xs,ys,xm,ym;
118: TSGetApplicationContext(ts,&user);
119: TSGetTime(ts,&t);
120: TSGetSolution(ts,&X);
122: DMGetCoordinateDM(user->da,&cda);
123: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
124: DMGetCoordinates(user->da,&gc);
125: DMDAVecGetArrayRead(cda,gc,&coors);
126: DMDAVecGetArrayRead(user->da,X,&p);
127: for (i=xs; i < xs+xm; i++) {
128: for (j=ys; j < ys+ym; j++) {
129: if (coors[j][i].y < 5) sum += p[j][i];
130: }
131: }
132: DMDAVecRestoreArrayRead(cda,gc,&coors);
133: DMDAVecRestoreArrayRead(user->da,X,&p);
134: MPI_Allreduce(&sum,&asum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ts));
135: PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %f = %f\n",(double)t,(double)(asum));
136: if (sum < 1.0e-2) {
137: TSSetConvergedReason(ts,TS_CONVERGED_USER);
138: PetscPrintf(PETSC_COMM_WORLD,"Exiting TS as the integral of PDF is almost zero\n");
139: }
140: return(0);
141: }
143: PetscErrorCode ini_bou(Vec X,AppCtx* user)
144: {
146: DM cda;
147: DMDACoor2d **coors;
148: PetscScalar **p;
149: Vec gc;
150: PetscInt i,j;
151: PetscInt xs,ys,xm,ym,M,N;
152: PetscScalar xi,yi;
153: PetscScalar sigmax=user->sigmax,sigmay=user->sigmay;
154: PetscScalar rho =user->rho;
155: PetscScalar muy=user->muy,mux;
156: PetscMPIInt rank;
157: PetscScalar sum;
160: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
161: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
162: user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
163: DMGetCoordinateDM(user->da,&cda);
164: DMGetCoordinates(user->da,&gc);
165: DMDAVecGetArray(cda,gc,&coors);
166: DMDAVecGetArray(user->da,X,&p);
167: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
169: /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable
170: muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point
171: in the y-direction. We only modify mux here
172: */
173: mux = user->mux = coors[0][M/2+10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */
174: if (user->nonoiseinitial) {
175: for (i=xs; i < xs+xm; i++) {
176: for (j=ys; j < ys+ym; j++) {
177: xi = coors[j][i].x; yi = coors[j][i].y;
178: if ((xi == mux) && (yi == muy)) {
179: p[j][i] = 1.0;
180: }
181: }
182: }
183: } else {
184: /* Change PM_min accordingly */
185: user->PM_min = user->Pmax*PetscSinScalar(mux);
186: for (i=xs; i < xs+xm; i++) {
187: for (j=ys; j < ys+ym; j++) {
188: xi = coors[j][i].x; yi = coors[j][i].y;
189: 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)));
190: }
191: }
192: }
193: DMDAVecRestoreArray(cda,gc,&coors);
194: DMDAVecRestoreArray(user->da,X,&p);
195: VecSum(X,&sum);
196: VecScale(X,1.0/sum);
197: return(0);
198: }
200: /* First advection term */
201: PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user)
202: {
203: PetscScalar f,fpos,fneg;
205: f = (y - user->ws);
206: fpos = PetscMax(f,0);
207: fneg = PetscMin(f,0);
208: if (user->st_width == 1) {
209: *p1 = fpos*(p[j][i] - p[j][i-1])/user->dx + fneg*(p[j][i+1] - p[j][i])/user->dx;
210: } else if (user->st_width == 2) {
211: *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);
212: } else if (user->st_width == 3) {
213: *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);
214: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils");
215: return(0);
216: }
218: /* Second advection term */
219: PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user)
220: {
221: PetscScalar f,fpos,fneg;
223: f = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x));
224: fpos = PetscMax(f,0);
225: fneg = PetscMin(f,0);
226: if (user->st_width == 1) {
227: *p2 = fpos*(p[j][i] - p[j-1][i])/user->dy + fneg*(p[j+1][i] - p[j][i])/user->dy;
228: } else if (user->st_width ==2) {
229: *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);
230: } else if (user->st_width == 3) {
231: *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);
232: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils");
233: return(0);
234: }
236: /* Diffusion term */
237: PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user)
238: {
240: if (user->st_width == 1) {
241: *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy));
242: } else if (user->st_width == 2) {
243: *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));
244: } else if (user->st_width == 3) {
245: *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));
246: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils");
247: return(0);
248: }
250: PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
251: {
253: AppCtx *user=(AppCtx*)ctx;
254: DM cda;
255: DMDACoor2d **coors;
256: PetscScalar **p,**f,**pdot;
257: PetscInt i,j;
258: PetscInt xs,ys,xm,ym,M,N;
259: Vec localX,gc,localXdot;
260: PetscScalar p_adv1,p_adv2,p_diff;
263: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
264: DMGetCoordinateDM(user->da,&cda);
265: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
267: DMGetLocalVector(user->da,&localX);
268: DMGetLocalVector(user->da,&localXdot);
270: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
271: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
272: DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);
273: DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);
275: DMGetCoordinatesLocal(user->da,&gc);
277: DMDAVecGetArrayRead(cda,gc,&coors);
278: DMDAVecGetArrayRead(user->da,localX,&p);
279: DMDAVecGetArrayRead(user->da,localXdot,&pdot);
280: DMDAVecGetArray(user->da,F,&f);
282: user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda));
283: for (i=xs; i < xs+xm; i++) {
284: for (j=ys; j < ys+ym; j++) {
285: adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);
286: adv2(p,coors[j][i].x,i,j,N,&p_adv2,user);
287: diffuse(p,i,j,t,&p_diff,user);
288: f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
289: }
290: }
291: DMDAVecRestoreArrayRead(user->da,localX,&p);
292: DMDAVecRestoreArrayRead(user->da,localX,&pdot);
293: DMRestoreLocalVector(user->da,&localX);
294: DMRestoreLocalVector(user->da,&localXdot);
295: DMDAVecRestoreArray(user->da,F,&f);
296: DMDAVecRestoreArrayRead(cda,gc,&coors);
298: return(0);
299: }
301: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx)
302: {
304: AppCtx *user=(AppCtx*)ctx;
305: DM cda;
306: DMDACoor2d **coors;
307: PetscInt i,j;
308: PetscInt xs,ys,xm,ym,M,N;
309: Vec gc;
310: PetscScalar val[5],xi,yi;
311: MatStencil row,col[5];
312: PetscScalar c1,c3,c5,c1pos,c1neg,c3pos,c3neg;
315: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
316: DMGetCoordinateDM(user->da,&cda);
317: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
319: DMGetCoordinatesLocal(user->da,&gc);
320: DMDAVecGetArrayRead(cda,gc,&coors);
321: for (i=xs; i < xs+xm; i++) {
322: for (j=ys; j < ys+ym; j++) {
323: PetscInt nc = 0;
324: xi = coors[j][i].x; yi = coors[j][i].y;
325: row.i = i; row.j = j;
326: c1 = (yi-user->ws)/user->dx;
327: c1pos = PetscMax(c1,0);
328: c1neg = PetscMin(c1,0);
329: c3 = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi))/user->dy;
330: c3pos = PetscMax(c3,0);
331: c3neg = PetscMin(c3,0);
332: c5 = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy);
333: col[nc].i = i-1; col[nc].j = j; val[nc++] = c1pos;
334: col[nc].i = i+1; col[nc].j = j; val[nc++] = -c1neg;
335: col[nc].i = i; col[nc].j = j-1; val[nc++] = c3pos + c5;
336: col[nc].i = i; col[nc].j = j+1; val[nc++] = -c3neg + c5;
337: col[nc].i = i; col[nc].j = j; val[nc++] = -c1pos + c1neg -c3pos + c3neg -2*c5 -a;
338: MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);
339: }
340: }
341: DMDAVecRestoreArrayRead(cda,gc,&coors);
343: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
344: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
345: if (J != Jpre) {
346: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
347: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
348: }
349: return(0);
350: }
354: PetscErrorCode Parameter_settings(AppCtx *user)
355: {
357: PetscBool flg;
361: /* Set default parameters */
362: user->ws = 1.0;
363: user->H = 5.0;
364: user->Pmax = 2.1;
365: user->PM_min = 1.0;
366: user->lambda = 0.1;
367: user->q = 1.0;
368: user->mux = PetscAsinScalar(user->PM_min/user->Pmax);
369: user->sigmax = 0.1;
370: user->sigmay = 0.1;
371: user->rho = 0.0;
372: user->xmin = -PETSC_PI;
373: user->xmax = PETSC_PI;
374: user->bx = DM_BOUNDARY_PERIODIC;
375: user->by = DM_BOUNDARY_MIRROR;
376: user->nonoiseinitial = PETSC_FALSE;
378: /*
379: ymin of -3 seems to let the unstable solution move up and leave a zero in its wake
380: with an ymin of -1 the wake is never exactly zero
381: */
382: user->ymin = -3.0;
383: user->ymax = 10.0;
384: user->st_width = 1;
386: PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg);
387: PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg);
388: PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg);
389: PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg);
390: PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg);
391: PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg);
392: PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg);
393: PetscOptionsGetScalar(NULL,NULL,"-sigmax",&user->sigmax,&flg);
394: PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg);
395: if (flg == 0) {
396: user->muy = user->ws;
397: }
398: PetscOptionsGetScalar(NULL,NULL,"-sigmay",&user->sigmay,&flg);
399: PetscOptionsGetScalar(NULL,NULL,"-rho",&user->rho,&flg);
400: PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg);
401: PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg);
402: PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg);
403: PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg);
404: PetscOptionsGetInt(NULL,NULL,"-stencil_width",&user->st_width,&flg);
405: PetscOptionsGetEnum(NULL,NULL,"-bx",DMBoundaryTypes,(PetscEnum*)&user->bx,&flg);
406: PetscOptionsGetEnum(NULL,NULL,"-by",DMBoundaryTypes,(PetscEnum*)&user->by,&flg);
407: PetscOptionsGetBool(NULL,NULL,"-nonoiseinitial",&user->nonoiseinitial,&flg);
409: return(0);
410: }
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*/