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;
61: PetscInitialize(&argc,&argv,"petscopt_ex7", help);
63: /* Get physics and time parameters */
64: Parameter_settings(&user);
65: /* Create a 2D DA with dof = 1 */
66: 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);
67: DMSetFromOptions(user.da);
68: DMSetUp(user.da);
69: /* Set x and y coordinates */
70: DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0,0);
71: DMDASetCoordinateName(user.da,0,"X - the angle");
72: DMDASetCoordinateName(user.da,1,"Y - the speed");
74: /* Get global vector x from DM */
75: DMCreateGlobalVector(user.da,&x);
77: ini_bou(x,&user);
78: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer);
79: VecView(x,viewer);
80: PetscViewerDestroy(&viewer);
82: TSCreate(PETSC_COMM_WORLD,&ts);
83: TSSetDM(ts,user.da);
84: TSSetProblemType(ts,TS_NONLINEAR);
85: TSSetType(ts,TSARKIMEX);
86: TSSetIFunction(ts,NULL,IFunction,&user);
87: /* TSSetIJacobian(ts,NULL,NULL,IJacobian,&user); */
88: TSSetApplicationContext(ts,&user);
89: TSSetTimeStep(ts,.005);
90: TSSetFromOptions(ts);
91: TSSetPostStep(ts,PostStep);
92: TSSolve(ts,x);
94: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fin_x",FILE_MODE_WRITE,&viewer);
95: VecView(x,viewer);
96: PetscViewerDestroy(&viewer);
98: VecDestroy(&x);
99: DMDestroy(&user.da);
100: TSDestroy(&ts);
101: PetscFinalize();
102: return 0;
103: }
105: PetscErrorCode PostStep(TS ts)
106: {
107: Vec X,gc;
108: AppCtx *user;
109: PetscScalar sum = 0,asum;
110: PetscReal t,**p;
111: DMDACoor2d **coors;
112: DM cda;
113: PetscInt i,j,xs,ys,xm,ym;
115: TSGetApplicationContext(ts,&user);
116: TSGetTime(ts,&t);
117: TSGetSolution(ts,&X);
119: DMGetCoordinateDM(user->da,&cda);
120: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
121: DMGetCoordinates(user->da,&gc);
122: DMDAVecGetArrayRead(cda,gc,&coors);
123: DMDAVecGetArrayRead(user->da,X,&p);
124: for (i=xs; i < xs+xm; i++) {
125: for (j=ys; j < ys+ym; j++) {
126: if (coors[j][i].y < 5) sum += p[j][i];
127: }
128: }
129: DMDAVecRestoreArrayRead(cda,gc,&coors);
130: DMDAVecRestoreArrayRead(user->da,X,&p);
131: MPI_Allreduce(&sum,&asum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)ts));
132: PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %f = %f\n",(double)t,(double)(asum));
133: if (sum < 1.0e-2) {
134: TSSetConvergedReason(ts,TS_CONVERGED_USER);
135: PetscPrintf(PETSC_COMM_WORLD,"Exiting TS as the integral of PDF is almost zero\n");
136: }
137: return 0;
138: }
140: PetscErrorCode ini_bou(Vec X,AppCtx* user)
141: {
142: DM cda;
143: DMDACoor2d **coors;
144: PetscScalar **p;
145: Vec gc;
146: PetscInt i,j;
147: PetscInt xs,ys,xm,ym,M,N;
148: PetscScalar xi,yi;
149: PetscScalar sigmax=user->sigmax,sigmay=user->sigmay;
150: PetscScalar rho =user->rho;
151: PetscScalar muy=user->muy,mux;
152: PetscMPIInt rank;
153: PetscScalar sum;
156: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
157: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
158: user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
159: DMGetCoordinateDM(user->da,&cda);
160: DMGetCoordinates(user->da,&gc);
161: DMDAVecGetArray(cda,gc,&coors);
162: DMDAVecGetArray(user->da,X,&p);
163: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
165: /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable
166: muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point
167: in the y-direction. We only modify mux here
168: */
169: mux = user->mux = coors[0][M/2+10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */
170: if (user->nonoiseinitial) {
171: for (i=xs; i < xs+xm; i++) {
172: for (j=ys; j < ys+ym; j++) {
173: xi = coors[j][i].x; yi = coors[j][i].y;
174: if ((xi == mux) && (yi == muy)) {
175: p[j][i] = 1.0;
176: }
177: }
178: }
179: } else {
180: /* Change PM_min accordingly */
181: user->PM_min = user->Pmax*PetscSinScalar(mux);
182: for (i=xs; i < xs+xm; i++) {
183: for (j=ys; j < ys+ym; j++) {
184: xi = coors[j][i].x; yi = coors[j][i].y;
185: 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)));
186: }
187: }
188: }
189: DMDAVecRestoreArray(cda,gc,&coors);
190: DMDAVecRestoreArray(user->da,X,&p);
191: VecSum(X,&sum);
192: VecScale(X,1.0/sum);
193: return 0;
194: }
196: /* First advection term */
197: PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user)
198: {
199: PetscScalar f,fpos,fneg;
200: f = (y - user->ws);
201: fpos = PetscMax(f,0);
202: fneg = PetscMin(f,0);
203: if (user->st_width == 1) {
204: *p1 = fpos*(p[j][i] - p[j][i-1])/user->dx + fneg*(p[j][i+1] - p[j][i])/user->dx;
205: } else if (user->st_width == 2) {
206: *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);
207: } else if (user->st_width == 3) {
208: *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);
209: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils");
210: return 0;
211: }
213: /* Second advection term */
214: PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user)
215: {
216: PetscScalar f,fpos,fneg;
217: f = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x));
218: fpos = PetscMax(f,0);
219: fneg = PetscMin(f,0);
220: if (user->st_width == 1) {
221: *p2 = fpos*(p[j][i] - p[j-1][i])/user->dy + fneg*(p[j+1][i] - p[j][i])/user->dy;
222: } else if (user->st_width ==2) {
223: *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);
224: } else if (user->st_width == 3) {
225: *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);
226: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils");
227: return 0;
228: }
230: /* Diffusion term */
231: PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user)
232: {
234: if (user->st_width == 1) {
235: *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy));
236: } else if (user->st_width == 2) {
237: *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));
238: } else if (user->st_width == 3) {
239: *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));
240: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"No support for wider stencils");
241: return 0;
242: }
244: PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
245: {
246: AppCtx *user=(AppCtx*)ctx;
247: DM cda;
248: DMDACoor2d **coors;
249: PetscScalar **p,**f,**pdot;
250: PetscInt i,j;
251: PetscInt xs,ys,xm,ym,M,N;
252: Vec localX,gc,localXdot;
253: PetscScalar p_adv1,p_adv2,p_diff;
256: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
257: DMGetCoordinateDM(user->da,&cda);
258: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
260: DMGetLocalVector(user->da,&localX);
261: DMGetLocalVector(user->da,&localXdot);
263: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
264: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
265: DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);
266: DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);
268: DMGetCoordinatesLocal(user->da,&gc);
270: DMDAVecGetArrayRead(cda,gc,&coors);
271: DMDAVecGetArrayRead(user->da,localX,&p);
272: DMDAVecGetArrayRead(user->da,localXdot,&pdot);
273: DMDAVecGetArray(user->da,F,&f);
275: user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda));
276: for (i=xs; i < xs+xm; i++) {
277: for (j=ys; j < ys+ym; j++) {
278: adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);
279: adv2(p,coors[j][i].x,i,j,N,&p_adv2,user);
280: diffuse(p,i,j,t,&p_diff,user);
281: f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
282: }
283: }
284: DMDAVecRestoreArrayRead(user->da,localX,&p);
285: DMDAVecRestoreArrayRead(user->da,localX,&pdot);
286: DMRestoreLocalVector(user->da,&localX);
287: DMRestoreLocalVector(user->da,&localXdot);
288: DMDAVecRestoreArray(user->da,F,&f);
289: DMDAVecRestoreArrayRead(cda,gc,&coors);
291: return 0;
292: }
294: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx)
295: {
296: AppCtx *user=(AppCtx*)ctx;
297: DM cda;
298: DMDACoor2d **coors;
299: PetscInt i,j;
300: PetscInt xs,ys,xm,ym,M,N;
301: Vec gc;
302: PetscScalar val[5],xi,yi;
303: MatStencil row,col[5];
304: PetscScalar c1,c3,c5,c1pos,c1neg,c3pos,c3neg;
307: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
308: DMGetCoordinateDM(user->da,&cda);
309: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
311: DMGetCoordinatesLocal(user->da,&gc);
312: DMDAVecGetArrayRead(cda,gc,&coors);
313: for (i=xs; i < xs+xm; i++) {
314: for (j=ys; j < ys+ym; j++) {
315: PetscInt nc = 0;
316: xi = coors[j][i].x; yi = coors[j][i].y;
317: row.i = i; row.j = j;
318: c1 = (yi-user->ws)/user->dx;
319: c1pos = PetscMax(c1,0);
320: c1neg = PetscMin(c1,0);
321: c3 = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi))/user->dy;
322: c3pos = PetscMax(c3,0);
323: c3neg = PetscMin(c3,0);
324: c5 = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy);
325: col[nc].i = i-1; col[nc].j = j; val[nc++] = c1pos;
326: col[nc].i = i+1; col[nc].j = j; val[nc++] = -c1neg;
327: col[nc].i = i; col[nc].j = j-1; val[nc++] = c3pos + c5;
328: col[nc].i = i; col[nc].j = j+1; val[nc++] = -c3neg + c5;
329: col[nc].i = i; col[nc].j = j; val[nc++] = -c1pos + c1neg -c3pos + c3neg -2*c5 -a;
330: MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);
331: }
332: }
333: DMDAVecRestoreArrayRead(cda,gc,&coors);
335: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
336: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
337: if (J != Jpre) {
338: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
339: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
340: }
341: return 0;
342: }
344: PetscErrorCode Parameter_settings(AppCtx *user)
345: {
346: PetscBool flg;
350: /* Set default parameters */
351: user->ws = 1.0;
352: user->H = 5.0;
353: user->Pmax = 2.1;
354: user->PM_min = 1.0;
355: user->lambda = 0.1;
356: user->q = 1.0;
357: user->mux = PetscAsinScalar(user->PM_min/user->Pmax);
358: user->sigmax = 0.1;
359: user->sigmay = 0.1;
360: user->rho = 0.0;
361: user->xmin = -PETSC_PI;
362: user->xmax = PETSC_PI;
363: user->bx = DM_BOUNDARY_PERIODIC;
364: user->by = DM_BOUNDARY_MIRROR;
365: user->nonoiseinitial = PETSC_FALSE;
367: /*
368: ymin of -3 seems to let the unstable solution move up and leave a zero in its wake
369: with an ymin of -1 the wake is never exactly zero
370: */
371: user->ymin = -3.0;
372: user->ymax = 10.0;
373: user->st_width = 1;
375: PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg);
376: PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg);
377: PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg);
378: PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg);
379: PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg);
380: PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg);
381: PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg);
382: PetscOptionsGetScalar(NULL,NULL,"-sigmax",&user->sigmax,&flg);
383: PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg);
384: if (flg == 0) {
385: user->muy = user->ws;
386: }
387: PetscOptionsGetScalar(NULL,NULL,"-sigmay",&user->sigmay,&flg);
388: PetscOptionsGetScalar(NULL,NULL,"-rho",&user->rho,&flg);
389: PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg);
390: PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg);
391: PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg);
392: PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg);
393: PetscOptionsGetInt(NULL,NULL,"-stencil_width",&user->st_width,&flg);
394: PetscOptionsGetEnum(NULL,NULL,"-bx",DMBoundaryTypes,(PetscEnum*)&user->bx,&flg);
395: PetscOptionsGetEnum(NULL,NULL,"-by",DMBoundaryTypes,(PetscEnum*)&user->by,&flg);
396: PetscOptionsGetBool(NULL,NULL,"-nonoiseinitial",&user->nonoiseinitial,&flg);
398: return 0;
399: }
401: /*TEST
403: build:
404: requires: !complex !single
406: test:
407: args: -ts_max_steps 2
408: localrunfiles: petscopt_ex7
410: test:
411: suffix: 2
412: args: -ts_max_steps 2 -snes_mf_operator
413: output_file: output/ex7_1.out
414: localrunfiles: petscopt_ex7
415: timeoutfactor: 2
417: test:
418: suffix: 3
419: args: -ts_max_steps 2 -snes_mf -pc_type none
420: output_file: output/ex7_1.out
421: localrunfiles: petscopt_ex7
422: timeoutfactor: 2
424: TEST*/