Actual source code: ex6.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;
  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: {
 53:   Vec            x;  /* Solution vector */
 54:   TS             ts;   /* Time-stepping context */
 55:   AppCtx         user; /* Application context */
 56:   Mat            J;
 57:   PetscViewer    viewer;

 59:   PetscInitialize(&argc,&argv,"petscopt_ex6", help);
 60:   /* Get physics and time parameters */
 61:   Parameter_settings(&user);
 62:   /* Create a 2D DA with dof = 1 */
 63:   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);
 64:   DMSetFromOptions(user.da);
 65:   DMSetUp(user.da);
 66:   /* Set x and y coordinates */
 67:   DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0.0,1.0);

 69:   /* Get global vector x from DM  */
 70:   DMCreateGlobalVector(user.da,&x);

 72:   ini_bou(x,&user);
 73:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer);
 74:   VecView(x,viewer);
 75:   PetscViewerDestroy(&viewer);

 77:   /* Get Jacobian matrix structure from the da */
 78:   DMSetMatType(user.da,MATAIJ);
 79:   DMCreateMatrix(user.da,&J);

 81:   TSCreate(PETSC_COMM_WORLD,&ts);
 82:   TSSetProblemType(ts,TS_NONLINEAR);
 83:   TSSetIFunction(ts,NULL,IFunction,&user);
 84:   TSSetIJacobian(ts,J,J,IJacobian,&user);
 85:   TSSetApplicationContext(ts,&user);
 86:   TSSetMaxTime(ts,user.tmax);
 87:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
 88:   TSSetTime(ts,user.t0);
 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:   MatDestroy(&J);
100:   DMDestroy(&user.da);
101:   TSDestroy(&ts);
102:   PetscFinalize();
103:   return 0;
104: }

106: PetscErrorCode PostStep(TS ts)
107: {
108:   Vec            X;
109:   AppCtx         *user;
110:   PetscScalar    sum;
111:   PetscReal      t;

113:   TSGetApplicationContext(ts,&user);
114:   TSGetTime(ts,&t);
115:   TSGetSolution(ts,&X);
116:   VecSum(X,&sum);
117:   PetscPrintf(PETSC_COMM_WORLD,"sum(p)*dw*dtheta at t = %3.2f = %3.6f\n",(double)t,(double)sum*user->dx*user->dy);
118:   return 0;
119: }

121: PetscErrorCode ini_bou(Vec X,AppCtx* user)
122: {
123:   DM             cda;
124:   DMDACoor2d     **coors;
125:   PetscScalar    **p;
126:   Vec            gc;
127:   PetscInt       i,j;
128:   PetscInt       xs,ys,xm,ym,M,N;
129:   PetscScalar    xi,yi;
130:   PetscScalar    sigmax=user->sigmax,sigmay=user->sigmay;
131:   PetscScalar    rho   =user->rho;
132:   PetscScalar    mux   =user->mux,muy=user->muy;
133:   PetscMPIInt    rank;

136:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
137:   DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
138:   user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
139:   DMGetCoordinateDM(user->da,&cda);
140:   DMGetCoordinates(user->da,&gc);
141:   DMDAVecGetArray(cda,gc,&coors);
142:   DMDAVecGetArray(user->da,X,&p);
143:   DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);
144:   for (i=xs; i < xs+xm; i++) {
145:     for (j=ys; j < ys+ym; j++) {
146:       xi = coors[j][i].x; yi = coors[j][i].y;
147:       if (i == 0 || j == 0 || i == M-1 || j == N-1) p[j][i] = 0.0;
148:       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)));
149:     }
150:   }
151:   /*  p[N/2+N%2][M/2+M%2] = 1/(user->dx*user->dy); */

153:   DMDAVecRestoreArray(cda,gc,&coors);
154:   DMDAVecRestoreArray(user->da,X,&p);
155:   return 0;
156: }

158: /* First advection term */
159: PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user)
160: {
161:   PetscScalar f;
162:   /*  PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */
163:   /*  if (i > 2 && i < M-2) {
164:     v1 = (y-user->ws)*(p[j][i-2] - p[j][i-3])/user->dx;
165:     v2 = (y-user->ws)*(p[j][i-1] - p[j][i-2])/user->dx;
166:     v3 = (y-user->ws)*(p[j][i] - p[j][i-1])/user->dx;
167:     v4 = (y-user->ws)*(p[j][i+1] - p[j][i])/user->dx;
168:     v5 = (y-user->ws)*(p[j][i+1] - p[j][i+2])/user->dx;

170:     s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
171:     s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
172:     s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;

174:     *p1 = 0.1*s1 + 0.6*s2 + 0.3*s3;
175:     } else *p1 = 0.0; */
176:   f   =  (y - user->ws);
177:   *p1 = f*(p[j][i+1] - p[j][i-1])/(2*user->dx);
178:   return 0;
179: }

181: /* Second advection term */
182: PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user)
183: {
184:   PetscScalar f;
185:   /*  PetscScalar v1,v2,v3,v4,v5,s1,s2,s3; */
186:   /*  if (j > 2 && j < N-2) {
187:     v1 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-2][i] - p[j-3][i])/user->dy;
188:     v2 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j-1][i] - p[j-2][i])/user->dy;
189:     v3 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j][i] - p[j-1][i])/user->dy;
190:     v4 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+1][i] - p[j][i])/user->dy;
191:     v5 = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*sin(x))*(p[j+2][i] - p[j+1][i])/user->dy;

193:     s1 = v1/3.0 - (7.0/6.0)*v2 + (11.0/6.0)*v3;
194:     s2 =-v2/6.0 + (5.0/6.0)*v3 + (1.0/3.0)*v4;
195:     s3 = v3/3.0 + (5.0/6.0)*v4 - (1.0/6.0)*v5;

197:     *p2 = 0.1*s1 + 0.6*s2 + 0.3*s3;
198:     } else *p2 = 0.0; */
199:   f   = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x));
200:   *p2 = f*(p[j+1][i] - p[j-1][i])/(2*user->dy);
201:   return 0;
202: }

204: /* Diffusion term */
205: PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user)
206: {

209:   *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy));
210:   return 0;
211: }

213: PetscErrorCode BoundaryConditions(PetscScalar **p,DMDACoor2d **coors,PetscInt i,PetscInt j,PetscInt M, PetscInt N,PetscScalar **f,AppCtx *user)
214: {
215:   PetscScalar fwc,fthetac;
216:   PetscScalar w=coors[j][i].y,theta=coors[j][i].x;

219:   if (user->bc == 0) { /* Natural boundary condition */
220:     f[j][i] = p[j][i];
221:   } else { /* Steady state boundary condition */
222:     fthetac = user->ws/(2*user->H)*(user->PM_min - user->Pmax*PetscSinScalar(theta));
223:     fwc = (w*w/2.0 - user->ws*w);
224:     if (i == 0 && j == 0) { /* left bottom corner */
225:       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;
226:     } else if (i == 0 && j == N-1) { /* right bottom corner */
227:       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;
228:     } else if (i == M-1 && j == 0) { /* left top corner */
229:       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;
230:     } else if (i == M-1 && j == N-1) { /* right top corner */
231:       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;
232:     } else if (i == 0) { /* Bottom edge */
233:       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);
234:     } else if (i == M-1) { /* Top edge */
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-1][i])/(2*user->dy);
236:     } else if (j == 0) { /* Left edge */
237:       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);
238:     } else if (j == N-1) { /* Right edge */
239:       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);
240:     }
241:   }
242:   return 0;
243: }

245: PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
246: {
247:   AppCtx         *user=(AppCtx*)ctx;
248:   DM             cda;
249:   DMDACoor2d     **coors;
250:   PetscScalar    **p,**f,**pdot;
251:   PetscInt       i,j;
252:   PetscInt       xs,ys,xm,ym,M,N;
253:   Vec            localX,gc,localXdot;
254:   PetscScalar    p_adv1,p_adv2,p_diff;

257:   DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
258:   DMGetCoordinateDM(user->da,&cda);
259:   DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);

261:   DMGetLocalVector(user->da,&localX);
262:   DMGetLocalVector(user->da,&localXdot);

264:   DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);
265:   DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);
266:   DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);
267:   DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);

269:   DMGetCoordinatesLocal(user->da,&gc);

271:   DMDAVecGetArrayRead(cda,gc,&coors);
272:   DMDAVecGetArrayRead(user->da,localX,&p);
273:   DMDAVecGetArrayRead(user->da,localXdot,&pdot);
274:   DMDAVecGetArray(user->da,F,&f);

276:   user->disper_coe = PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda));
277:   for (i=xs; i < xs+xm; i++) {
278:     for (j=ys; j < ys+ym; j++) {
279:       if (i == 0 || j == 0 || i == M-1 || j == N-1) {
280:         BoundaryConditions(p,coors,i,j,M,N,f,user);
281:       } else {
282:         adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);
283:         adv2(p,coors[j][i].x,i,j,N,&p_adv2,user);
284:         diffuse(p,i,j,t,&p_diff,user);
285:         f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
286:       }
287:     }
288:   }
289:   DMDAVecRestoreArrayRead(user->da,localX,&p);
290:   DMDAVecRestoreArrayRead(user->da,localX,&pdot);
291:   DMRestoreLocalVector(user->da,&localX);
292:   DMRestoreLocalVector(user->da,&localXdot);
293:   DMDAVecRestoreArray(user->da,F,&f);
294:   DMDAVecRestoreArrayRead(cda,gc,&coors);

296:   return 0;
297: }

299: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx)
300: {
301:   AppCtx         *user=(AppCtx*)ctx;
302:   DM             cda;
303:   DMDACoor2d     **coors;
304:   PetscInt       i,j;
305:   PetscInt       xs,ys,xm,ym,M,N;
306:   Vec            gc;
307:   PetscScalar    val[5],xi,yi;
308:   MatStencil     row,col[5];
309:   PetscScalar    c1,c3,c5;

312:   DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
313:   DMGetCoordinateDM(user->da,&cda);
314:   DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);

316:   DMGetCoordinatesLocal(user->da,&gc);
317:   DMDAVecGetArrayRead(cda,gc,&coors);
318:   for (i=xs; i < xs+xm; i++) {
319:     for (j=ys; j < ys+ym; j++) {
320:       PetscInt nc = 0;
321:       xi = coors[j][i].x; yi = coors[j][i].y;
322:       row.i = i; row.j = j;
323:       if (i == 0 || j == 0 || i == M-1 || j == N-1) {
324:         if (user->bc == 0) {
325:           col[nc].i = i; col[nc].j = j; val[nc++] = 1.0;
326:         } else {
327:           PetscScalar fthetac,fwc;
328:           fthetac = user->ws/(2*user->H)*(user->PM_min - user->Pmax*PetscSinScalar(xi));
329:           fwc     = (yi*yi/2.0 - user->ws*yi);
330:           if (i==0 && j==0) {
331:             col[nc].i = i+1; col[nc].j = j;   val[nc++] = fwc/user->dx;
332:             col[nc].i = i;   col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy;
333:             col[nc].i = i;   col[nc].j = j;   val[nc++] = -fwc/user->dx + fthetac + user->disper_coe/user->dy;
334:           } else if (i==0 && j == N-1) {
335:             col[nc].i = i+1; col[nc].j = j;   val[nc++] = fwc/user->dx;
336:             col[nc].i = i;   col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy;
337:             col[nc].i = i;   col[nc].j = j;   val[nc++] = -fwc/user->dx + fthetac - user->disper_coe/user->dy;
338:           } else if (i== M-1 && 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 == M-1 && 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==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/(2*user->dy);
349:             col[nc].i = i;   col[nc].j = j-1; val[nc++] =  user->disper_coe/(2*user->dy);
350:             col[nc].i = i;   col[nc].j = j;   val[nc++] = -fwc/user->dx + fthetac;
351:           } else if (i == M-1) {
352:             col[nc].i = i-1; col[nc].j = j;   val[nc++] = -fwc/user->dx;
353:             col[nc].i = i;   col[nc].j = j+1; val[nc++] = -user->disper_coe/(2*user->dy);
354:             col[nc].i = i;   col[nc].j = j-1; val[nc++] =  user->disper_coe/(2*user->dy);
355:             col[nc].i = i;   col[nc].j = j;   val[nc++] = fwc/user->dx + fthetac;
356:           } else if (j==0) {
357:             col[nc].i = i+1; col[nc].j = j;   val[nc++] = fwc/(2*user->dx);
358:             col[nc].i = i-1; col[nc].j = j;   val[nc++] = -fwc/(2*user->dx);
359:             col[nc].i = i;   col[nc].j = j+1; val[nc++] = -user->disper_coe/user->dy;
360:             col[nc].i = i;   col[nc].j = j;   val[nc++] = user->disper_coe/user->dy + fthetac;
361:           } else if (j == N-1) {
362:             col[nc].i = i+1; col[nc].j = j;   val[nc++] = fwc/(2*user->dx);
363:             col[nc].i = i-1; col[nc].j = j;   val[nc++] = -fwc/(2*user->dx);
364:             col[nc].i = i;   col[nc].j = j-1; val[nc++] = user->disper_coe/user->dy;
365:             col[nc].i = i;   col[nc].j = j;   val[nc++] = -user->disper_coe/user->dy + fthetac;
366:           }
367:         }
368:       } else {
369:         c1        = (yi-user->ws)/(2*user->dx);
370:         c3        = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi))/(2*user->dy);
371:         c5        = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy);
372:         col[nc].i = i-1; col[nc].j = j;   val[nc++] = c1;
373:         col[nc].i = i+1; col[nc].j = j;   val[nc++] = -c1;
374:         col[nc].i = i;   col[nc].j = j-1; val[nc++] = c3 + c5;
375:         col[nc].i = i;   col[nc].j = j+1; val[nc++] = -c3 + c5;
376:         col[nc].i = i;   col[nc].j = j;   val[nc++] = -2*c5 -a;
377:       }
378:       MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);
379:     }
380:   }
381:   DMDAVecRestoreArrayRead(cda,gc,&coors);

383:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
384:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
385:   if (J != Jpre) {
386:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
387:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
388:   }
389:   return 0;
390: }

392: PetscErrorCode Parameter_settings(AppCtx *user)
393: {
394:   PetscBool      flg;


398:   /* Set default parameters */
399:   user->ws     = 1.0;
400:   user->H      = 5.0;  user->Pmax   = 2.1;
401:   user->PM_min = 1.0;  user->lambda = 0.1;
402:   user->q      = 1.0;  user->mux    = PetscAsinScalar(user->PM_min/user->Pmax);
403:   user->sigmax = 0.1;
404:   user->sigmay = 0.1;  user->rho  = 0.0;
405:   user->t0     = 0.0;  user->tmax = 2.0;
406:   user->xmin   = -1.0; user->xmax = 10.0;
407:   user->ymin   = -1.0; user->ymax = 10.0;
408:   user->bc     = 0;

410:   PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg);
411:   PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg);
412:   PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg);
413:   PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg);
414:   PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg);
415:   PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg);
416:   PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg);
417:   PetscOptionsGetScalar(NULL,NULL,"-sigmax",&user->sigmax,&flg);
418:   PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg);
419:   PetscOptionsGetScalar(NULL,NULL,"-sigmay",&user->sigmay,&flg);
420:   PetscOptionsGetScalar(NULL,NULL,"-rho",&user->rho,&flg);
421:   PetscOptionsGetScalar(NULL,NULL,"-t0",&user->t0,&flg);
422:   PetscOptionsGetScalar(NULL,NULL,"-tmax",&user->tmax,&flg);
423:   PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg);
424:   PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg);
425:   PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg);
426:   PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg);
427:   PetscOptionsGetInt(NULL,NULL,"-bc_type",&user->bc,&flg);
428:   user->muy = user->ws;
429:   return 0;
430: }

432: /*TEST

434:    build:
435:       requires: !complex

437:    test:
438:       args: -nox -ts_max_steps 2
439:       localrunfiles: petscopt_ex6
440:       timeoutfactor: 4

442: TEST*/