Actual source code: ex6.c

petsc-3.12.0 2019-09-29
Report Typos and Errors
  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*/