xref: /petsc/src/ts/tutorials/power_grid/ex7.c (revision 0baf8eba40dbc839082666f9f7396a225d6f663c)
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 Neumman using mirror values
7 
8    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.
9    x_t = (y - ws)  y_t = (ws/2H)*(Pm - Pmax*sin(x))
10 
11 */
12 
13 #include <petscdm.h>
14 #include <petscdmda.h>
15 #include <petscts.h>
16 
17 /*
18    User-defined data structures and routines
19 */
20 typedef struct {
21   PetscScalar    ws;         /* Synchronous speed */
22   PetscScalar    H;          /* Inertia constant */
23   PetscScalar    D;          /* Damping constant */
24   PetscScalar    Pmax;       /* Maximum power output of generator */
25   PetscScalar    PM_min;     /* Mean mechanical power input */
26   PetscScalar    lambda;     /* correlation time */
27   PetscScalar    q;          /* noise strength */
28   PetscScalar    mux;        /* Initial average angle */
29   PetscScalar    sigmax;     /* Standard deviation of initial angle */
30   PetscScalar    muy;        /* Average speed */
31   PetscScalar    sigmay;     /* standard deviation of initial speed */
32   PetscScalar    rho;        /* Cross-correlation coefficient */
33   PetscScalar    xmin;       /* left boundary of angle */
34   PetscScalar    xmax;       /* right boundary of angle */
35   PetscScalar    ymin;       /* bottom boundary of speed */
36   PetscScalar    ymax;       /* top boundary of speed */
37   PetscScalar    dx;         /* x step size */
38   PetscScalar    dy;         /* y step size */
39   PetscScalar    disper_coe; /* Dispersion coefficient */
40   DM             da;
41   PetscInt       st_width; /* Stencil width */
42   DMBoundaryType bx;       /* x boundary type */
43   DMBoundaryType by;       /* y boundary type */
44   PetscBool      nonoiseinitial;
45 } AppCtx;
46 
47 PetscErrorCode Parameter_settings(AppCtx *);
48 PetscErrorCode ini_bou(Vec, AppCtx *);
49 PetscErrorCode IFunction(TS, PetscReal, Vec, Vec, Vec, void *);
50 PetscErrorCode IJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
51 PetscErrorCode PostStep(TS);
52 
53 int main(int argc, char **argv)
54 {
55   Vec         x;    /* Solution vector */
56   TS          ts;   /* Time-stepping context */
57   AppCtx      user; /* Application context */
58   PetscViewer viewer;
59 
60   PetscFunctionBeginUser;
61   PetscCall(PetscInitialize(&argc, &argv, "petscopt_ex7", help));
62 
63   /* Get physics and time parameters */
64   PetscCall(Parameter_settings(&user));
65   /* Create a 2D DA with dof = 1 */
66   PetscCall(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   PetscCall(DMSetFromOptions(user.da));
68   PetscCall(DMSetUp(user.da));
69   /* Set x and y coordinates */
70   PetscCall(DMDASetUniformCoordinates(user.da, user.xmin, user.xmax, user.ymin, user.ymax, 0, 0));
71   PetscCall(DMDASetCoordinateName(user.da, 0, "X - the angle"));
72   PetscCall(DMDASetCoordinateName(user.da, 1, "Y - the speed"));
73 
74   /* Get global vector x from DM  */
75   PetscCall(DMCreateGlobalVector(user.da, &x));
76 
77   PetscCall(ini_bou(x, &user));
78   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "ini_x", FILE_MODE_WRITE, &viewer));
79   PetscCall(VecView(x, viewer));
80   PetscCall(PetscViewerDestroy(&viewer));
81 
82   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
83   PetscCall(TSSetDM(ts, user.da));
84   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
85   PetscCall(TSSetType(ts, TSARKIMEX));
86   PetscCall(TSSetIFunction(ts, NULL, IFunction, &user));
87   /*  PetscCall(TSSetIJacobian(ts,NULL,NULL,IJacobian,&user));  */
88   PetscCall(TSSetApplicationContext(ts, &user));
89   PetscCall(TSSetTimeStep(ts, .005));
90   PetscCall(TSSetFromOptions(ts));
91   PetscCall(TSSetPostStep(ts, PostStep));
92   PetscCall(TSSolve(ts, x));
93 
94   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "fin_x", FILE_MODE_WRITE, &viewer));
95   PetscCall(VecView(x, viewer));
96   PetscCall(PetscViewerDestroy(&viewer));
97 
98   PetscCall(VecDestroy(&x));
99   PetscCall(DMDestroy(&user.da));
100   PetscCall(TSDestroy(&ts));
101   PetscCall(PetscFinalize());
102   return 0;
103 }
104 
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;
114 
115   PetscFunctionBegin;
116   PetscCall(TSGetApplicationContext(ts, &user));
117   PetscCall(TSGetTime(ts, &t));
118   PetscCall(TSGetSolution(ts, &X));
119 
120   PetscCall(DMGetCoordinateDM(user->da, &cda));
121   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
122   PetscCall(DMGetCoordinates(user->da, &gc));
123   PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
124   PetscCall(DMDAVecGetArrayRead(user->da, X, &p));
125   for (i = xs; i < xs + xm; i++) {
126     for (j = ys; j < ys + ym; j++) {
127       if (coors[j][i].y < 5) sum += p[j][i];
128     }
129   }
130   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
131   PetscCall(DMDAVecRestoreArrayRead(user->da, X, &p));
132   PetscCallMPI(MPIU_Allreduce(&sum, &asum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
133   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "sum(p)*dw*dtheta at t = %f = %f\n", (double)t, (double)asum));
134   if (sum < 1.0e-2) {
135     PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_USER));
136     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Exiting TS as the integral of PDF is almost zero\n"));
137   }
138   PetscFunctionReturn(PETSC_SUCCESS);
139 }
140 
141 PetscErrorCode ini_bou(Vec X, AppCtx *user)
142 {
143   DM            cda;
144   DMDACoor2d  **coors;
145   PetscScalar **p;
146   Vec           gc;
147   PetscInt      i, j;
148   PetscInt      xs, ys, xm, ym, M, N;
149   PetscScalar   xi, yi;
150   PetscScalar   sigmax = user->sigmax, sigmay = user->sigmay;
151   PetscScalar   rho = user->rho;
152   PetscScalar   muy = user->muy, mux;
153   PetscMPIInt   rank;
154   PetscScalar   sum;
155 
156   PetscFunctionBeginUser;
157   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
158   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
159   user->dx = (user->xmax - user->xmin) / (M - 1);
160   user->dy = (user->ymax - user->ymin) / (N - 1);
161   PetscCall(DMGetCoordinateDM(user->da, &cda));
162   PetscCall(DMGetCoordinates(user->da, &gc));
163   PetscCall(DMDAVecGetArray(cda, gc, &coors));
164   PetscCall(DMDAVecGetArray(user->da, X, &p));
165   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
166 
167   /* mux and muy need to be grid points in the x and y-direction otherwise the solution goes unstable
168      muy is set by choosing the y domain, no. of grid points along y-direction so that muy is a grid point
169      in the y-direction. We only modify mux here
170   */
171   mux = user->mux = coors[0][M / 2 + 10].x; /* For -pi < x < pi, this should be some angle between 0 and pi/2 */
172   if (user->nonoiseinitial) {
173     for (i = xs; i < xs + xm; i++) {
174       for (j = ys; j < ys + ym; j++) {
175         xi = coors[j][i].x;
176         yi = coors[j][i].y;
177         if ((xi == mux) && (yi == muy)) p[j][i] = 1.0;
178       }
179     }
180   } else {
181     /* Change PM_min accordingly */
182     user->PM_min = user->Pmax * PetscSinScalar(mux);
183     for (i = xs; i < xs + xm; i++) {
184       for (j = ys; j < ys + ym; j++) {
185         xi = coors[j][i].x;
186         yi = coors[j][i].y;
187         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)));
188       }
189     }
190   }
191   PetscCall(DMDAVecRestoreArray(cda, gc, &coors));
192   PetscCall(DMDAVecRestoreArray(user->da, X, &p));
193   PetscCall(VecSum(X, &sum));
194   PetscCall(VecScale(X, 1.0 / sum));
195   PetscFunctionReturn(PETSC_SUCCESS);
196 }
197 
198 /* First advection term */
199 PetscErrorCode adv1(PetscScalar **p, PetscScalar y, PetscInt i, PetscInt j, PetscInt M, PetscScalar *p1, AppCtx *user)
200 {
201   PetscScalar f, fpos, fneg;
202 
203   PetscFunctionBegin;
204   f    = (y - user->ws);
205   fpos = PetscMax(f, 0);
206   fneg = PetscMin(f, 0);
207   if (user->st_width == 1) {
208     *p1 = fpos * (p[j][i] - p[j][i - 1]) / user->dx + fneg * (p[j][i + 1] - p[j][i]) / user->dx;
209   } else if (user->st_width == 2) {
210     *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);
211   } else if (user->st_width == 3) {
212     *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);
213   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
214   PetscFunctionReturn(PETSC_SUCCESS);
215 }
216 
217 /* Second advection term */
218 PetscErrorCode adv2(PetscScalar **p, PetscScalar x, PetscInt i, PetscInt j, PetscInt N, PetscScalar *p2, AppCtx *user)
219 {
220   PetscScalar f, fpos, fneg;
221 
222   PetscFunctionBegin;
223   f    = (user->ws / (2 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(x));
224   fpos = PetscMax(f, 0);
225   fneg = PetscMin(f, 0);
226   if (user->st_width == 1) {
227     *p2 = fpos * (p[j][i] - p[j - 1][i]) / user->dy + fneg * (p[j + 1][i] - p[j][i]) / user->dy;
228   } else if (user->st_width == 2) {
229     *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);
230   } else if (user->st_width == 3) {
231     *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);
232   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
233   PetscFunctionReturn(PETSC_SUCCESS);
234 }
235 
236 /* Diffusion term */
237 PetscErrorCode diffuse(PetscScalar **p, PetscInt i, PetscInt j, PetscReal t, PetscScalar *p_diff, AppCtx *user)
238 {
239   PetscFunctionBeginUser;
240   if (user->st_width == 1) {
241     *p_diff = user->disper_coe * ((p[j - 1][i] - 2 * p[j][i] + p[j + 1][i]) / (user->dy * user->dy));
242   } else if (user->st_width == 2) {
243     *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));
244   } else if (user->st_width == 3) {
245     *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));
246   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No support for wider stencils");
247   PetscFunctionReturn(PETSC_SUCCESS);
248 }
249 
250 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
251 {
252   AppCtx       *user = (AppCtx *)ctx;
253   DM            cda;
254   DMDACoor2d  **coors;
255   PetscScalar **p, **f, **pdot;
256   PetscInt      i, j;
257   PetscInt      xs, ys, xm, ym, M, N;
258   Vec           localX, gc, localXdot;
259   PetscScalar   p_adv1 = 0.0, p_adv2 = 0.0, p_diff; /* initialize to prevent incorrect compiler warnings */
260 
261   PetscFunctionBeginUser;
262   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
263   PetscCall(DMGetCoordinateDM(user->da, &cda));
264   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
265 
266   PetscCall(DMGetLocalVector(user->da, &localX));
267   PetscCall(DMGetLocalVector(user->da, &localXdot));
268 
269   PetscCall(DMGlobalToLocalBegin(user->da, X, INSERT_VALUES, localX));
270   PetscCall(DMGlobalToLocalEnd(user->da, X, INSERT_VALUES, localX));
271   PetscCall(DMGlobalToLocalBegin(user->da, Xdot, INSERT_VALUES, localXdot));
272   PetscCall(DMGlobalToLocalEnd(user->da, Xdot, INSERT_VALUES, localXdot));
273 
274   PetscCall(DMGetCoordinatesLocal(user->da, &gc));
275 
276   PetscCall(DMDAVecGetArrayRead(cda, gc, &coors));
277   PetscCall(DMDAVecGetArrayRead(user->da, localX, &p));
278   PetscCall(DMDAVecGetArrayRead(user->da, localXdot, &pdot));
279   PetscCall(DMDAVecGetArray(user->da, F, &f));
280 
281   user->disper_coe = PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda));
282   for (i = xs; i < xs + xm; i++) {
283     for (j = ys; j < ys + ym; j++) {
284       PetscCall(adv1(p, coors[j][i].y, i, j, M, &p_adv1, user));
285       PetscCall(adv2(p, coors[j][i].x, i, j, N, &p_adv2, user));
286       PetscCall(diffuse(p, i, j, t, &p_diff, user));
287       f[j][i] = -p_adv1 - p_adv2 + p_diff - pdot[j][i];
288     }
289   }
290   PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &p));
291   PetscCall(DMDAVecRestoreArrayRead(user->da, localX, &pdot));
292   PetscCall(DMRestoreLocalVector(user->da, &localX));
293   PetscCall(DMRestoreLocalVector(user->da, &localXdot));
294   PetscCall(DMDAVecRestoreArray(user->da, F, &f));
295   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
296   PetscFunctionReturn(PETSC_SUCCESS);
297 }
298 
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, c1pos, c1neg, c3pos, c3neg;
310 
311   PetscFunctionBeginUser;
312   PetscCall(DMDAGetInfo(user->da, NULL, &M, &N, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
313   PetscCall(DMGetCoordinateDM(user->da, &cda));
314   PetscCall(DMDAGetCorners(cda, &xs, &ys, 0, &xm, &ym, 0));
315 
316   PetscCall(DMGetCoordinatesLocal(user->da, &gc));
317   PetscCall(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;
322       yi          = coors[j][i].y;
323       row.i       = i;
324       row.j       = j;
325       c1          = (yi - user->ws) / user->dx;
326       c1pos       = PetscMax(c1, 0);
327       c1neg       = PetscMin(c1, 0);
328       c3          = (user->ws / (2.0 * user->H)) * (user->PM_min - user->Pmax * PetscSinScalar(xi)) / user->dy;
329       c3pos       = PetscMax(c3, 0);
330       c3neg       = PetscMin(c3, 0);
331       c5          = (PetscPowScalar((user->lambda * user->ws) / (2 * user->H), 2) * user->q * (1.0 - PetscExpScalar(-t / user->lambda))) / (user->dy * user->dy);
332       col[nc].i   = i - 1;
333       col[nc].j   = j;
334       val[nc++]   = c1pos;
335       col[nc].i   = i + 1;
336       col[nc].j   = j;
337       val[nc++]   = -c1neg;
338       col[nc].i   = i;
339       col[nc].j   = j - 1;
340       val[nc++]   = c3pos + c5;
341       col[nc].i   = i;
342       col[nc].j   = j + 1;
343       val[nc++]   = -c3neg + c5;
344       col[nc].i   = i;
345       col[nc].j   = j;
346       val[nc++]   = -c1pos + c1neg - c3pos + c3neg - 2 * c5 - a;
347       PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
348     }
349   }
350   PetscCall(DMDAVecRestoreArrayRead(cda, gc, &coors));
351 
352   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
353   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
354   if (J != Jpre) {
355     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
356     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
357   }
358   PetscFunctionReturn(PETSC_SUCCESS);
359 }
360 
361 PetscErrorCode Parameter_settings(AppCtx *user)
362 {
363   PetscBool flg;
364 
365   PetscFunctionBeginUser;
366   /* Set default parameters */
367   user->ws             = 1.0;
368   user->H              = 5.0;
369   user->Pmax           = 2.1;
370   user->PM_min         = 1.0;
371   user->lambda         = 0.1;
372   user->q              = 1.0;
373   user->mux            = PetscAsinScalar(user->PM_min / user->Pmax);
374   user->sigmax         = 0.1;
375   user->sigmay         = 0.1;
376   user->rho            = 0.0;
377   user->xmin           = -PETSC_PI;
378   user->xmax           = PETSC_PI;
379   user->bx             = DM_BOUNDARY_PERIODIC;
380   user->by             = DM_BOUNDARY_MIRROR;
381   user->nonoiseinitial = PETSC_FALSE;
382 
383   /*
384      ymin of -3 seems to let the unstable solution move up and leave a zero in its wake
385      with an ymin of -1 the wake is never exactly zero
386   */
387   user->ymin     = -3.0;
388   user->ymax     = 10.0;
389   user->st_width = 1;
390 
391   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ws", &user->ws, &flg));
392   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Inertia", &user->H, &flg));
393   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-Pmax", &user->Pmax, &flg));
394   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-PM_min", &user->PM_min, &flg));
395   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-lambda", &user->lambda, &flg));
396   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-q", &user->q, &flg));
397   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mux", &user->mux, &flg));
398   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmax", &user->sigmax, &flg));
399   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-muy", &user->muy, &flg));
400   if (flg == 0) user->muy = user->ws;
401   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-sigmay", &user->sigmay, &flg));
402   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-rho", &user->rho, &flg));
403   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmin", &user->xmin, &flg));
404   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-xmax", &user->xmax, &flg));
405   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymin", &user->ymin, &flg));
406   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-ymax", &user->ymax, &flg));
407   PetscCall(PetscOptionsGetInt(NULL, NULL, "-stencil_width", &user->st_width, &flg));
408   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-bx", DMBoundaryTypes, (PetscEnum *)&user->bx, &flg));
409   PetscCall(PetscOptionsGetEnum(NULL, NULL, "-by", DMBoundaryTypes, (PetscEnum *)&user->by, &flg));
410   PetscCall(PetscOptionsGetBool(NULL, NULL, "-nonoiseinitial", &user->nonoiseinitial, &flg));
411   PetscFunctionReturn(PETSC_SUCCESS);
412 }
413 
414 /*TEST
415 
416    build:
417       requires: !complex !single
418 
419    test:
420       args: -ts_max_steps 2
421       localrunfiles: petscopt_ex7
422 
423    test:
424       suffix: 2
425       args: -ts_max_steps 2 -snes_mf_operator
426       output_file: output/ex7_1.out
427       localrunfiles: petscopt_ex7
428       timeoutfactor: 2
429 
430    test:
431       suffix: 3
432       args: -ts_max_steps 2 -snes_mf -pc_type none
433       output_file: output/ex7_1.out
434       localrunfiles: petscopt_ex7
435       timeoutfactor: 2
436 
437 TEST*/
438