xref: /petsc/src/ts/tutorials/power_grid/ex8.c (revision 4e278199b78715991f5c71ebbd945c1489263e6c)
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 
6    Boundary conditions:
7    Zero dirichlet in y using ghosted values
8    Periodic in x
9 
10    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.
11    x_t = (y - ws)
12    y_t = (ws/2H)*(Pm - Pmax*sin(x) - D*(w - ws))
13 
14    In this example, we can see the effect of a fault, that zeroes the electrical power output
15    Pmax*sin(x), on the PDF. The fault on/off times can be controlled by options -tf and -tcl respectively.
16 
17 */
18 
19 #include <petscdm.h>
20 #include <petscdmda.h>
21 #include <petscts.h>
22 
23 /*
24    User-defined data structures and routines
25 */
26 typedef struct {
27   PetscScalar ws;   /* Synchronous speed */
28   PetscScalar H;    /* Inertia constant */
29   PetscScalar D;    /* Damping constant */
30   PetscScalar Pmax,Pmax_s; /* Maximum power output of generator */
31   PetscScalar PM_min; /* Mean mechanical power input */
32   PetscScalar lambda; /* correlation time */
33   PetscScalar q;      /* noise strength */
34   PetscScalar mux;    /* Initial average angle */
35   PetscScalar sigmax; /* Standard deviation of initial angle */
36   PetscScalar muy;    /* Average speed */
37   PetscScalar sigmay; /* standard deviation of initial speed */
38   PetscScalar rho;    /* Cross-correlation coefficient */
39   PetscScalar xmin;   /* left boundary of angle */
40   PetscScalar xmax;   /* right boundary of angle */
41   PetscScalar ymin;   /* bottom boundary of speed */
42   PetscScalar ymax;   /* top boundary of speed */
43   PetscScalar dx;     /* x step size */
44   PetscScalar dy;     /* y step size */
45   PetscScalar disper_coe; /* Dispersion coefficient */
46   DM          da;
47   PetscInt    st_width; /* Stencil width */
48   DMBoundaryType bx; /* x boundary type */
49   DMBoundaryType by; /* y boundary type */
50   PetscReal        tf,tcl; /* Fault incidence and clearing times */
51 } AppCtx;
52 
53 PetscErrorCode Parameter_settings(AppCtx*);
54 PetscErrorCode ini_bou(Vec,AppCtx*);
55 PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
56 PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
57 PetscErrorCode PostStep(TS);
58 
59 int main(int argc, char **argv)
60 {
61   PetscErrorCode ierr;
62   Vec            x;  /* Solution vector */
63   TS             ts;   /* Time-stepping context */
64   AppCtx         user; /* Application context */
65   PetscViewer    viewer;
66 
67   ierr = PetscInitialize(&argc,&argv,"petscopt_ex8", help);if (ierr) return ierr;
68 
69   /* Get physics and time parameters */
70   ierr = Parameter_settings(&user);CHKERRQ(ierr);
71   /* Create a 2D DA with dof = 1 */
72   ierr = 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);CHKERRQ(ierr);
73   ierr = DMSetFromOptions(user.da);CHKERRQ(ierr);
74   ierr = DMSetUp(user.da);CHKERRQ(ierr);
75   /* Set x and y coordinates */
76   ierr = DMDASetUniformCoordinates(user.da,user.xmin,user.xmax,user.ymin,user.ymax,0,0);CHKERRQ(ierr);
77   ierr = DMDASetCoordinateName(user.da,0,"X - the angle");CHKERRQ(ierr);
78   ierr = DMDASetCoordinateName(user.da,1,"Y - the speed");CHKERRQ(ierr);
79 
80   /* Get global vector x from DM  */
81   ierr = DMCreateGlobalVector(user.da,&x);CHKERRQ(ierr);
82 
83   ierr = ini_bou(x,&user);CHKERRQ(ierr);
84   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"ini_x",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
85   ierr = VecView(x,viewer);CHKERRQ(ierr);
86   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
87 
88   ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
89   ierr = TSSetDM(ts,user.da);CHKERRQ(ierr);
90   ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
91   ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);
92   ierr = TSSetIFunction(ts,NULL,IFunction,&user);CHKERRQ(ierr);
93   /*  ierr = TSSetIJacobian(ts,NULL,NULL,IJacobian,&user);CHKERRQ(ierr); */
94   ierr = TSSetApplicationContext(ts,&user);CHKERRQ(ierr);
95   ierr = TSSetTimeStep(ts,.005);CHKERRQ(ierr);
96   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
97   ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr);
98   ierr = TSSolve(ts,x);CHKERRQ(ierr);
99 
100   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"fin_x",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
101   ierr = VecView(x,viewer);CHKERRQ(ierr);
102   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
103 
104   ierr = VecDestroy(&x);CHKERRQ(ierr);
105   ierr = DMDestroy(&user.da);CHKERRQ(ierr);
106   ierr = TSDestroy(&ts);CHKERRQ(ierr);
107   ierr = PetscFinalize();
108   return ierr;
109 }
110 
111 PetscErrorCode PostStep(TS ts)
112 {
113   PetscErrorCode ierr;
114   Vec            X;
115   AppCtx         *user;
116   PetscReal      t;
117   PetscScalar    asum;
118 
119   PetscFunctionBegin;
120   ierr = TSGetApplicationContext(ts,&user);CHKERRQ(ierr);
121   ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
122   ierr = TSGetSolution(ts,&X);CHKERRQ(ierr);
123   /*
124   if (t >= .2) {
125     ierr = TSGetSolution(ts,&X);CHKERRQ(ierr);
126     ierr = VecView(X,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr);
127     exit(0);
128      results in initial conditions after fault in binaryoutput
129   }*/
130 
131   if ((t > user->tf) && (t < user->tcl)) user->Pmax = 0.0; /* A short-circuit that drives the electrical power output (Pmax*sin(delta)) to zero */
132   else user->Pmax = user->Pmax_s;
133 
134   ierr = VecSum(X,&asum);CHKERRQ(ierr);
135   ierr = PetscPrintf(PETSC_COMM_WORLD,"sum(p) at t = %f = %f\n",(double)t,(double)(asum));CHKERRQ(ierr);
136   PetscFunctionReturn(0);
137 }
138 
139 PetscErrorCode ini_bou(Vec X,AppCtx* user)
140 {
141   PetscErrorCode ierr;
142   DM             cda;
143   DMDACoor2d     **coors;
144   PetscScalar    **p;
145   Vec            gc;
146   PetscInt       M,N,Ir,J;
147   PetscMPIInt    rank;
148 
149   PetscFunctionBeginUser;
150   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
151   ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
152   user->dx = (user->xmax - user->xmin)/(M-1); user->dy = (user->ymax - user->ymin)/(N-1);
153   ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);
154   ierr = DMGetCoordinates(user->da,&gc);CHKERRQ(ierr);
155   ierr = DMDAVecGetArrayRead(cda,gc,&coors);CHKERRQ(ierr);
156   ierr = DMDAVecGetArray(user->da,X,&p);CHKERRQ(ierr);
157 
158   /* Point mass at (mux,muy) */
159   ierr = PetscPrintf(PETSC_COMM_WORLD,"Original user->mux = %f, user->muy = %f\n",user->mux,user->muy);CHKERRQ(ierr);
160   ierr = DMDAGetLogicalCoordinate(user->da,user->mux,user->muy,0.0,&Ir,&J,NULL,&user->mux,&user->muy,NULL);CHKERRQ(ierr);
161   user->PM_min = user->Pmax*PetscSinScalar(user->mux);
162   ierr = PetscPrintf(PETSC_COMM_WORLD,"Corrected user->mux = %f, user->muy = %f user->PM_min = %f,user->dx = %f\n",user->mux,user->muy,user->PM_min,user->dx);CHKERRQ(ierr);
163   if (Ir > -1 && J > -1) {
164     p[J][Ir] = 1.0;
165   }
166 
167   ierr = DMDAVecRestoreArrayRead(cda,gc,&coors);CHKERRQ(ierr);
168   ierr = DMDAVecRestoreArray(user->da,X,&p);CHKERRQ(ierr);
169   PetscFunctionReturn(0);
170 }
171 
172 /* First advection term */
173 PetscErrorCode adv1(PetscScalar **p,PetscScalar y,PetscInt i,PetscInt j,PetscInt M,PetscScalar *p1,AppCtx *user)
174 {
175   PetscScalar f,fpos,fneg;
176   PetscFunctionBegin;
177   f   =  (y - user->ws);
178   fpos = PetscMax(f,0);
179   fneg = PetscMin(f,0);
180   if (user->st_width == 1) {
181     *p1 = fpos*(p[j][i] - p[j][i-1])/user->dx + fneg*(p[j][i+1] - p[j][i])/user->dx;
182   } else if (user->st_width == 2) {
183     *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);
184   } else if (user->st_width == 3) {
185     *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);
186   }
187   /* *p1 = f*(p[j][i+1] - p[j][i-1])/user->dx;*/
188   PetscFunctionReturn(0);
189 }
190 
191 /* Second advection term */
192 PetscErrorCode adv2(PetscScalar **p,PetscScalar x,PetscScalar y,PetscInt i,PetscInt j,PetscInt N,PetscScalar *p2,AppCtx *user)
193 {
194   PetscScalar f,fpos,fneg;
195   PetscFunctionBegin;
196   f   = (user->ws/(2*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(x) - user->D*(y - user->ws));
197   fpos = PetscMax(f,0);
198   fneg = PetscMin(f,0);
199   if (user->st_width == 1) {
200     *p2 = fpos*(p[j][i] - p[j-1][i])/user->dy + fneg*(p[j+1][i] - p[j][i])/user->dy;
201   } else if (user->st_width ==2) {
202     *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);
203   } else if (user->st_width == 3) {
204     *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);
205   }
206 
207   /* *p2 = f*(p[j+1][i] - p[j-1][i])/user->dy;*/
208   PetscFunctionReturn(0);
209 }
210 
211 /* Diffusion term */
212 PetscErrorCode diffuse(PetscScalar **p,PetscInt i,PetscInt j,PetscReal t,PetscScalar *p_diff,AppCtx * user)
213 {
214   PetscFunctionBeginUser;
215   if (user->st_width == 1) {
216     *p_diff = user->disper_coe*((p[j-1][i] - 2*p[j][i] + p[j+1][i])/(user->dy*user->dy));
217   } else if (user->st_width == 2) {
218     *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));
219   } else if (user->st_width == 3) {
220     *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));
221   }
222   PetscFunctionReturn(0);
223 }
224 
225 PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
226 {
227   PetscErrorCode ierr;
228   AppCtx         *user=(AppCtx*)ctx;
229   DM             cda;
230   DMDACoor2d     **coors;
231   PetscScalar    **p,**f,**pdot;
232   PetscInt       i,j;
233   PetscInt       xs,ys,xm,ym,M,N;
234   Vec            localX,gc,localXdot;
235   PetscScalar    p_adv1 = 0.0,p_adv2 = 0.0,p_diff = 0;
236   PetscScalar    diffuse1,gamma;
237 
238   PetscFunctionBeginUser;
239   ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
240   ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);
241   ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
242 
243   ierr = DMGetLocalVector(user->da,&localX);CHKERRQ(ierr);
244   ierr = DMGetLocalVector(user->da,&localXdot);CHKERRQ(ierr);
245 
246   ierr = DMGlobalToLocalBegin(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
247   ierr = DMGlobalToLocalEnd(user->da,X,INSERT_VALUES,localX);CHKERRQ(ierr);
248   ierr = DMGlobalToLocalBegin(user->da,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr);
249   ierr = DMGlobalToLocalEnd(user->da,Xdot,INSERT_VALUES,localXdot);CHKERRQ(ierr);
250 
251   ierr = DMGetCoordinatesLocal(user->da,&gc);CHKERRQ(ierr);
252 
253   ierr = DMDAVecGetArrayRead(cda,gc,&coors);CHKERRQ(ierr);
254   ierr = DMDAVecGetArrayRead(user->da,localX,&p);CHKERRQ(ierr);
255   ierr = DMDAVecGetArrayRead(user->da,localXdot,&pdot);CHKERRQ(ierr);
256   ierr = DMDAVecGetArray(user->da,F,&f);CHKERRQ(ierr);
257 
258   gamma = user->D*user->ws/(2*user->H);
259   diffuse1 = user->lambda*user->lambda*user->q/(user->lambda*gamma+1)*(1.0 - PetscExpScalar(-t*(gamma+1.0)/user->lambda));
260   user->disper_coe = user->ws*user->ws/(4*user->H*user->H)*diffuse1;
261 
262   for (i=xs; i < xs+xm; i++) {
263     for (j=ys; j < ys+ym; j++) {
264       ierr = adv1(p,coors[j][i].y,i,j,M,&p_adv1,user);CHKERRQ(ierr);
265       ierr = adv2(p,coors[j][i].x,coors[j][i].y,i,j,N,&p_adv2,user);CHKERRQ(ierr);
266       ierr = diffuse(p,i,j,t,&p_diff,user);CHKERRQ(ierr);
267       f[j][i] = -p_adv1 - p_adv2  + p_diff - pdot[j][i];
268     }
269   }
270   ierr = DMDAVecRestoreArrayRead(user->da,localX,&p);CHKERRQ(ierr);
271   ierr = DMDAVecRestoreArrayRead(user->da,localX,&pdot);CHKERRQ(ierr);
272   ierr = DMRestoreLocalVector(user->da,&localX);CHKERRQ(ierr);
273   ierr = DMRestoreLocalVector(user->da,&localXdot);CHKERRQ(ierr);
274   ierr = DMDAVecRestoreArray(user->da,F,&f);CHKERRQ(ierr);
275   ierr = DMDAVecRestoreArrayRead(cda,gc,&coors);CHKERRQ(ierr);
276 
277   PetscFunctionReturn(0);
278 }
279 
280 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ctx)
281 {
282   PetscErrorCode ierr;
283   AppCtx         *user=(AppCtx*)ctx;
284   DM             cda;
285   DMDACoor2d     **coors;
286   PetscInt       i,j;
287   PetscInt       xs,ys,xm,ym,M,N;
288   Vec            gc;
289   PetscScalar    val[5],xi,yi;
290   MatStencil     row,col[5];
291   PetscScalar    c1,c3,c5,c1pos,c1neg,c3pos,c3neg;
292 
293   PetscFunctionBeginUser;
294   ierr = DMDAGetInfo(user->da,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);CHKERRQ(ierr);
295   ierr = DMGetCoordinateDM(user->da,&cda);CHKERRQ(ierr);
296   ierr = DMDAGetCorners(cda,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr);
297 
298   ierr = DMGetCoordinatesLocal(user->da,&gc);CHKERRQ(ierr);
299   ierr = DMDAVecGetArrayRead(cda,gc,&coors);CHKERRQ(ierr);
300   for (i=xs; i < xs+xm; i++) {
301     for (j=ys; j < ys+ym; j++) {
302       PetscInt nc = 0;
303       xi = coors[j][i].x; yi = coors[j][i].y;
304       row.i = i; row.j = j;
305       c1        = (yi-user->ws)/user->dx;
306       c1pos    = PetscMax(c1,0);
307       c1neg    = PetscMin(c1,0);
308       c3        = (user->ws/(2.0*user->H))*(user->PM_min - user->Pmax*PetscSinScalar(xi) - user->D*(yi - user->ws))/user->dy;
309       c3pos    = PetscMax(c3,0);
310       c3neg    = PetscMin(c3,0);
311       c5        = (PetscPowScalar((user->lambda*user->ws)/(2*user->H),2)*user->q*(1.0-PetscExpScalar(-t/user->lambda)))/(user->dy*user->dy);
312       col[nc].i = i-1; col[nc].j = j;   val[nc++] = c1pos;
313       col[nc].i = i+1; col[nc].j = j;   val[nc++] = -c1neg;
314       col[nc].i = i;   col[nc].j = j-1; val[nc++] = c3pos + c5;
315       col[nc].i = i;   col[nc].j = j+1; val[nc++] = -c3neg + c5;
316       col[nc].i = i;   col[nc].j = j;   val[nc++] = -c1pos + c1neg -c3pos + c3neg -2*c5 -a;
317       ierr = MatSetValuesStencil(Jpre,1,&row,nc,col,val,INSERT_VALUES);CHKERRQ(ierr);
318     }
319   }
320   ierr = DMDAVecRestoreArrayRead(cda,gc,&coors);CHKERRQ(ierr);
321 
322   ierr =  MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
323   ierr = MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
324   if (J != Jpre) {
325     ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
326     ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
327   }
328   PetscFunctionReturn(0);
329 }
330 
331 PetscErrorCode Parameter_settings(AppCtx *user)
332 {
333   PetscErrorCode ierr;
334   PetscBool      flg;
335 
336   PetscFunctionBeginUser;
337 
338   /* Set default parameters */
339   user->ws     = 1.0;
340   user->H      = 5.0;
341   user->D      = 0.0;
342   user->Pmax = user->Pmax_s  = 2.1;
343   user->PM_min = 1.0;
344   user->lambda = 0.1;
345   user->q      = 1.0;
346   user->mux    = PetscAsinScalar(user->PM_min/user->Pmax);
347   user->sigmax = 0.1;
348   user->sigmay = 0.1;
349   user->rho    = 0.0;
350   user->xmin   = -PETSC_PI;
351   user->xmax   = PETSC_PI;
352   user->bx     = DM_BOUNDARY_PERIODIC;
353   user->by     = DM_BOUNDARY_GHOSTED;
354   user->tf = user->tcl = -1;
355   user->ymin   = -2.0;
356   user->ymax   = 2.0;
357   user->st_width = 1;
358 
359   ierr = PetscOptionsGetScalar(NULL,NULL,"-ws",&user->ws,&flg);CHKERRQ(ierr);
360   ierr = PetscOptionsGetScalar(NULL,NULL,"-Inertia",&user->H,&flg);CHKERRQ(ierr);
361   ierr = PetscOptionsGetScalar(NULL,NULL,"-D",&user->D,&flg);CHKERRQ(ierr);
362   ierr = PetscOptionsGetScalar(NULL,NULL,"-Pmax",&user->Pmax,&flg);CHKERRQ(ierr);
363   ierr = PetscOptionsGetScalar(NULL,NULL,"-PM_min",&user->PM_min,&flg);CHKERRQ(ierr);
364   ierr = PetscOptionsGetScalar(NULL,NULL,"-lambda",&user->lambda,&flg);CHKERRQ(ierr);
365   ierr = PetscOptionsGetScalar(NULL,NULL,"-q",&user->q,&flg);CHKERRQ(ierr);
366   ierr = PetscOptionsGetScalar(NULL,NULL,"-mux",&user->mux,&flg);CHKERRQ(ierr);
367   ierr = PetscOptionsGetScalar(NULL,NULL,"-muy",&user->muy,&flg);CHKERRQ(ierr);
368   if (flg == 0) {
369     user->muy = user->ws;
370   }
371   ierr = PetscOptionsGetScalar(NULL,NULL,"-xmin",&user->xmin,&flg);CHKERRQ(ierr);
372   ierr = PetscOptionsGetScalar(NULL,NULL,"-xmax",&user->xmax,&flg);CHKERRQ(ierr);
373   ierr = PetscOptionsGetScalar(NULL,NULL,"-ymin",&user->ymin,&flg);CHKERRQ(ierr);
374   ierr = PetscOptionsGetScalar(NULL,NULL,"-ymax",&user->ymax,&flg);CHKERRQ(ierr);
375   ierr = PetscOptionsGetInt(NULL,NULL,"-stencil_width",&user->st_width,&flg);CHKERRQ(ierr);
376   ierr = PetscOptionsGetEnum(NULL,NULL,"-bx",DMBoundaryTypes,(PetscEnum*)&user->bx,&flg);CHKERRQ(ierr);
377   ierr = PetscOptionsGetEnum(NULL,NULL,"-by",DMBoundaryTypes,(PetscEnum*)&user->by,&flg);CHKERRQ(ierr);
378   ierr = PetscOptionsGetReal(NULL,NULL,"-tf",&user->tf,&flg);CHKERRQ(ierr);
379   ierr = PetscOptionsGetReal(NULL,NULL,"-tcl",&user->tcl,&flg);CHKERRQ(ierr);
380   PetscFunctionReturn(0);
381 }
382 
383 /*TEST
384 
385    build:
386       requires: !complex x
387 
388    test:
389       args: -ts_max_steps 1
390       localrunfiles: petscopt_ex8
391       timeoutfactor: 3
392 
393 TEST*/
394