Actual source code: ex6.c
petsc-3.13.4 2020-08-01
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: x_t = (y - ws) y_t = (ws/2H)*(Pm - Pmax*sin(x))
7: Boundary conditions: -bc_type 0 => Zero dirichlet boundary
8: -bc_type 1 => Steady state boundary condition
9: Steady state boundary condition found by setting p_t = 0
10: */
12: #include <petscdm.h>
13: #include <petscdmda.h>
14: #include <petscts.h>
16: /*
17: User-defined data structures and routines
18: */
19: typedef struct {
20: PetscScalar ws; /* Synchronous speed */
21: PetscScalar H; /* Inertia constant */
22: PetscScalar D; /* Damping constant */
23: PetscScalar Pmax; /* Maximum power output of generator */
24: PetscScalar PM_min; /* Mean mechanical power input */
25: PetscScalar lambda; /* correlation time */
26: PetscScalar q; /* noise strength */
27: PetscScalar mux; /* Initial average angle */
28: PetscScalar sigmax; /* Standard deviation of initial angle */
29: PetscScalar muy; /* Average speed */
30: PetscScalar sigmay; /* standard deviation of initial speed */
31: PetscScalar rho; /* Cross-correlation coefficient */
32: PetscScalar t0; /* Initial time */
33: PetscScalar tmax; /* Final time */
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: PetscInt bc; /* Boundary conditions */
41: PetscScalar disper_coe; /* Dispersion coefficient */
42: DM da;
43: } AppCtx;
45: PetscErrorCode Parameter_settings(AppCtx*);
46: PetscErrorCode ini_bou(Vec,AppCtx*);
47: PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
48: PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
49: PetscErrorCode PostStep(TS);
51: int main(int argc, char **argv)
52: {
54: Vec x; /* Solution vector */
55: TS ts; /* Time-stepping context */
56: AppCtx user; /* Application context */
57: Mat J;
58: PetscViewer viewer;
60: PetscInitialize(&argc,&argv,"petscopt_ex6", help);if (ierr) return ierr;
61: /* Get physics and time parameters */
62: Parameter_settings(&user);
63: /* Create a 2D DA with dof = 1 */
64: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&user.da);
65: DMSetFromOptions(user.da);
66: DMSetUp(user.da);
67: /* Set x and y coordinates */
68: DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);
70: /* Get global vector x from DM */
71: DMCreateGlobalVector(user.da,&x);
73: ini_bou(x,&user);
74: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer);
75: VecView(x,viewer);
76: PetscViewerDestroy(&viewer);
78: /* Get Jacobian matrix structure from the da */
79: DMSetMatType(user.da,MATAIJ);
80: DMCreateMatrix(user.da,&J);
82: TSCreate(PETSC_COMM_WORLD,&ts);
83: TSSetProblemType(ts,TS_NONLINEAR);
84: TSSetIFunction(ts,NULL,IFunction,&user);
85: TSSetIJacobian(ts,J,J,IJacobian,&user);
86: TSSetApplicationContext(ts,&user);
87: TSSetMaxTime(ts,user.tmax);
88: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
89: TSSetTime(ts,user.t0);
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: MatDestroy(&J);
101: DMDestroy(&user.da);
102: TSDestroy(&ts);
103: PetscFinalize();
104: return ierr;
105: }
107: PetscErrorCode PostStep(TS ts)
108: {
110: Vec X;
111: AppCtx *user;
112: PetscScalar sum;
113: PetscReal t;
116: TSGetApplicationContext(ts,&user);
117: TSGetTime(ts,&t);
118: TSGetSolution(ts,&X);
119: VecSum(X,&sum);
120: PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %3.2f = %3.6f\n",(double)t,(double)sum*user->dx*user->dy);
121: return(0);
122: }
124: PetscErrorCode ini_bou(Vec X,AppCtx* user)
125: {
127: DM cda;
128: DMDACoor2d **coors;
129: PetscScalar **p;
130: Vec gc;
131: PetscInt i,j;
132: PetscInt xs,ys,xm,ym,M,N;
133: PetscScalar xi,yi;
134: PetscScalar sigmax=user->sigmax,sigmay=user->sigmay;
135: PetscScalar rho =user->rho;
136: PetscScalar mux =user->mux,muy=user->muy;
137: PetscMPIInt rank;
140: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
141: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
142: user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
143: DMGetCoordinateDM(user->da,&cda);
144: DMGetCoordinates(user->da,&gc);
145: DMDAVecGetArray(cda,gc,&coors);
146: DMDAVecGetArray(user->da,X,&p);
147: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
148: for (i=xs; i < xs+xm; i++) {
149: for (j=ys; j < ys+ym; j++) {
150: xi = coors[j][i].x; yi = coors[j][i].y;
151: if (i == 0 || j == 0 || i == M-1 || j == N-1) p[j][i] = 0.0;
152: else 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)));
153: }
154: }
155: /* p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */
157: DMDAVecRestoreArray(cda,gc,&coors);
158: DMDAVecRestoreArray(user->da,X,&p);
159: return(0);
160: }
162: /* First advection term */
163: PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user)
164: {
165: PetscScalar f;
166: /* PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */
168: /* if (i > 2 && i < M-2) {
169: v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx;
170: v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx;
171: v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx;
172: v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx;
173: v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx;
175: s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
176: s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
177: s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;
179: *p1 = 0.1*s1 + 0.6*s2 + 0.3*s3;
180: } else *p1 = 0.0; */
181: f = (y - user->ws);
182: *p1 = f*(p[j][i+1] - p[j][i-1])/(2*user->dx);
183: return(0);
184: }
186: /* Second advection term */
187: PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user)
188: {
189: PetscScalar f;
190: /* PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */
192: /* if (j > 2 && j < N-2) {
193: v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy;
194: v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy;
195: v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy;
196: v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy;
197: v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy;
199: s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
200: s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
201: s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;
203: *p2 = 0.1*s1 + 0.6*s2 + 0.3*s3;
204: } else *p2 = 0.0; */
205: f = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x));
206: *p2 = f*(p[j+1][i] - p[j-1][i])/(2*user->dy);
207: return(0);
208: }
210: /* Diffusion term */
211: PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user)
212: {
215: *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy));
216: return(0);
217: }
219: PetscErrorCode BoundaryConditions(PetscScalar **p,DMDACoor2d **coors,PetscInt i,PetscInt j,PetscInt M, PetscInt N,PetscScalar **f,AppCtx *user)
220: {
221: PetscScalar fwc,fthetac;
222: PetscScalar w=coors[j][i].y,theta=coors[j][i].x;
225: if (user->bc == 0) { /* Natural boundary condition */
226: f[j][i] = p[j][i];
227: } else { /* Steady state boundary condition */
228: fthetac = user->ws/(2*user->H)*(user->PM_min - user->Pmax*PetscSinScalar(theta));
229: fwc = (w*w/2.0 - user->ws*w);
230: if (i == 0 && j == 0) { /* left bottom corner */
231: f[j][i] = fwc*(p[j][i+1] - p[j][i])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j][i])/user->dy;
232: } else if (i == 0 && j == N-1) { /* right bottom corner */
233: f[j][i] = fwc*(p[j][i+1] - p[j][i])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j][i] - p[j-1][i])/user->dy;
234: } else if (i == M-1 && j == 0) { /* left top corner */
235: f[j][i] = fwc*(p[j][i] - p[j][i-1])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j][i])/user->dy;
236: } else if (i == M-1 && j == N-1) { /* right top corner */
237: f[j][i] = fwc*(p[j][i] - p[j][i-1])/user->dx + fthetac*p[j][i] - user->disper_coe*(p[j][i] - p[j-1][i])/user->dy;
238: } else if (i == 0) { /* Bottom edge */
239: f[j][i] = fwc*(p[j][i+1] - p[j][i])/(user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j-1][i])/(2*user->dy);
240: } else if (i == M-1) { /* Top edge */
241: f[j][i] = fwc*(p[j][i] - p[j][i-1])/(user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j-1][i])/(2*user->dy);
242: } else if (j == 0) { /* Left edge */
243: f[j][i] = fwc*(p[j][i+1] - p[j][i-1])/(2*user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j+1][i] - p[j][i])/(user->dy);
244: } else if (j == N-1) { /* Right edge */
245: f[j][i] = fwc*(p[j][i+1] - p[j][i-1])/(2*user->dx) + fthetac*p[j][i] - user->disper_coe*(p[j][i] - p[j-1][i])/(user->dy);
246: }
247: }
248: return(0);
249: }
251: PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
252: {
254: AppCtx *user=(AppCtx*)ctx;
255: DM cda;
256: DMDACoor2d **coors;
257: PetscScalar **p,**f,**pdot;
258: PetscInt i,j;
259: PetscInt xs,ys,xm,ym,M,N;
260: Vec localX,gc,localXdot;
261: PetscScalar p_adv1,p_adv2,p_diff;
264: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
265: DMGetCoordinateDM(user->da,&cda);
266: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
268: DMGetLocalVector(user->da,&localX);
269: DMGetLocalVector(user->da,&localXdot);
271: DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
272: DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
273: DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);
274: DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);
276: DMGetCoordinatesLocal(user->da,&gc);
278: DMDAVecGetArrayRead(cda,gc,&coors);
279: DMDAVecGetArrayRead(user->da,localX,&p);
280: DMDAVecGetArrayRead(user->da,localXdot,&pdot);
281: DMDAVecGetArray(user->da,F,&f);
283: user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda));
284: for (i=xs; i < xs+xm; i++) {
285: for (j=ys; j < ys+ym; j++) {
286: if (i == 0 || j == 0 || i == M-1 || j == N-1) {
287: BoundaryConditions(p,coors,i,j,M,N,f,user);
288: } else {
289: adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);
290: adv2(p,coors[j][i].x,i,j,N,&p_adv2,user);
291: diffuse(p,i,j,t,&p_diff,user);
292: f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
293: }
294: }
295: }
296: DMDAVecRestoreArrayRead(user->da,localX,&p);
297: DMDAVecRestoreArrayRead(user->da,localX,&pdot);
298: DMRestoreLocalVector(user->da,&localX);
299: DMRestoreLocalVector(user->da,&localXdot);
300: DMDAVecRestoreArray(user->da,F,&f);
301: DMDAVecRestoreArrayRead(cda,gc,&coors);
303: return(0);
304: }
306: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx)
307: {
309: AppCtx *user=(AppCtx*)ctx;
310: DM cda;
311: DMDACoor2d **coors;
312: PetscInt i,j;
313: PetscInt xs,ys,xm,ym,M,N;
314: Vec gc;
315: PetscScalar val[5],xi,yi;
316: MatStencil row,col[5];
317: PetscScalar c1,c3,c5;
320: DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
321: DMGetCoordinateDM(user->da,&cda);
322: DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
324: DMGetCoordinatesLocal(user->da,&gc);
325: DMDAVecGetArrayRead(cda,gc,&coors);
326: for (i=xs; i < xs+xm; i++) {
327: for (j=ys; j < ys+ym; j++) {
328: PetscInt nc = 0;
329: xi = coors[j][i].x; yi = coors[j][i].y;
330: row.i = i; row.j = j;
331: if (i == 0 || j == 0 || i == M-1 || j == N-1) {
332: if (user->bc == 0) {
333: col[nc].i = i; col[nc].j = j; val[nc++] = 1.0;
334: } else {
335: PetscScalar fthetac,fwc;
336: fthetac = user->ws/(2*user->H)*(user->PM_min - user->Pmax*PetscSinScalar(xi));
337: fwc = (yi*yi/2.0 - user->ws*yi);
338: if (i==0 && j==0) {
339: col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/user->dx;
340: col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy;
341: col[nc].i = i; col[nc].j = j; val[nc++] = -fwc/user->dx + fthetac + user->disper_coe/user->dy;
342: } else if (i==0 && j == N-1) {
343: col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/user->dx;
344: col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy;
345: col[nc].i = i; col[nc].j = j; val[nc++] = -fwc/user->dx + fthetac - user->disper_coe/user->dy;
346: } else if (i== M-1 && j == 0) {
347: col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/user->dx;
348: col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy;
349: col[nc].i = i; col[nc].j = j; val[nc++] = fwc/user->dx + fthetac + user->disper_coe/user->dy;
350: } else if (i == M-1 && j == N-1) {
351: col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/user->dx;
352: col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy;
353: col[nc].i = i; col[nc].j = j; val[nc++] = fwc/user->dx + fthetac - user->disper_coe/user->dy;
354: } else if (i==0) {
355: col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/user->dx;
356: col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/(2*user->dy);
357: col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/(2*user->dy);
358: col[nc].i = i; col[nc].j = j; val[nc++] = -fwc/user->dx + fthetac;
359: } else if (i == M-1) {
360: col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/user->dx;
361: col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/(2*user->dy);
362: col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/(2*user->dy);
363: col[nc].i = i; col[nc].j = j; val[nc++] = fwc/user->dx + fthetac;
364: } else if (j==0) {
365: col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/(2*user->dx);
366: col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/(2*user->dx);
367: col[nc].i = i; col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy;
368: col[nc].i = i; col[nc].j = j; val[nc++] = user->disper_coe/user->dy + fthetac;
369: } else if (j == N-1) {
370: col[nc].i = i+1; col[nc].j = j; val[nc++] = fwc/(2*user->dx);
371: col[nc].i = i-1; col[nc].j = j; val[nc++] = -fwc/(2*user->dx);
372: col[nc].i = i; col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy;
373: col[nc].i = i; col[nc].j = j; val[nc++] = -user->disper_coe/user->dy + fthetac;
374: }
375: }
376: } else {
377: c1 = (yi-user->ws)/(2*user->dx);
378: c3 = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi))/(2*user->dy);
379: c5 = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy);
380: col[nc].i = i-1; col[nc].j = j; val[nc++] = c1;
381: col[nc].i = i+1; col[nc].j = j; val[nc++] = -c1;
382: col[nc].i = i; col[nc].j = j-1; val[nc++] = c3 + c5;
383: col[nc].i = i; col[nc].j = j+1; val[nc++] = -c3 + c5;
384: col[nc].i = i; col[nc].j = j; val[nc++] = -2*c5 -a;
385: }
386: MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);
387: }
388: }
389: DMDAVecRestoreArrayRead(cda,gc,&coors);
391: MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
392: MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
393: if (J != Jpre) {
394: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
395: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
396: }
397: return(0);
398: }
400: PetscErrorCode Parameter_settings(AppCtx *user)
401: {
403: PetscBool flg;
407: /* Set default parameters */
408: user->ws = 1.0;
409: user->H = 5.0; user->Pmax = 2.1;
410: user->PM_min = 1.0; user->lambda = 0.1;
411: user->q = 1.0; user->mux = PetscAsinScalar(user->PM_min/user->Pmax);
412: user->sigmax = 0.1;
413: user->sigmay = 0.1; user->rho = 0.0;
414: user->t0 = 0.0; user->tmax = 2.0;
415: user->xmin = -1.0; user->xmax = 10.0;
416: user->ymin = -1.0; user->ymax = 10.0;
417: user->bc = 0;
418:
419: PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg);
420: PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg);
421: PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg);
422: PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg);
423: PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg);
424: PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg);
425: PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg);
426: PetscOptionsGetScalar(NULL,NULL,"-sigmax",&user->sigmax,&flg);
427: PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg);
428: PetscOptionsGetScalar(NULL,NULL,"-sigmay",&user->sigmay,&flg);
429: PetscOptionsGetScalar(NULL,NULL,"-rho",&user->rho,&flg);
430: PetscOptionsGetScalar(NULL,NULL,"-t0",&user->t0,&flg);
431: PetscOptionsGetScalar(NULL,NULL,"-tmax",&user->tmax,&flg);
432: PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg);
433: PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg);
434: PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg);
435: PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg);
436: PetscOptionsGetInt(NULL,NULL,"-bc_type",&user->bc,&flg);
437: user->muy = user->ws;
438: return(0);
439: }
442: /*TEST
444: build:
445: requires: !complex
447: test:
448: args: -nox -ts_max_steps 2
449: localrunfiles: petscopt_ex6
450: timeoutfactor: 4
452: TEST*/