xref: /petsc/src/ts/tutorials/power_grid/ex9opt.c (revision d5b43468fb8780a8feea140ccd6fa3e6a50411cc)
1 
2 static char help[] = "Basic equation for generator stability analysis.\n";
3 
4 /*F
5 
6 \begin{eqnarray}
7                  \frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
8                  \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
9 \end{eqnarray}
10 
11   Ensemble of initial conditions
12    ./ex2 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3      -ts_adapt_dt_max .01  -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
13 
14   Fault at .1 seconds
15    ./ex2           -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01  -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
16 
17   Initial conditions same as when fault is ended
18    ./ex2 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05  -ts_adapt_dt_max .01  -ts_monitor -ts_type rosw -pc_type lu -ksp_type preonly
19 
20 F*/
21 
22 /*
23    Include "petscts.h" so that we can use TS solvers.  Note that this
24    file automatically includes:
25      petscsys.h       - base PETSc routines   petscvec.h - vectors
26      petscmat.h - matrices
27      petscis.h     - index sets            petscksp.h - Krylov subspace methods
28      petscviewer.h - viewers               petscpc.h  - preconditioners
29      petscksp.h   - linear solvers
30 */
31 
32 #include <petsctao.h>
33 #include <petscts.h>
34 
35 typedef struct {
36   TS          ts;
37   PetscScalar H, D, omega_b, omega_s, Pmax, Pm, E, V, X, u_s, c;
38   PetscInt    beta;
39   PetscReal   tf, tcl, dt;
40 } AppCtx;
41 
42 PetscErrorCode FormFunction(Tao, Vec, PetscReal *, void *);
43 PetscErrorCode FormGradient(Tao, Vec, Vec, void *);
44 
45 /*
46      Defines the ODE passed to the ODE solver
47 */
48 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
49 {
50   PetscScalar       *f, Pmax;
51   const PetscScalar *u;
52 
53   PetscFunctionBegin;
54   /*  The next three lines allow us to access the entries of the vectors directly */
55   PetscCall(VecGetArrayRead(U, &u));
56   PetscCall(VecGetArray(F, &f));
57   if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
58   else Pmax = ctx->Pmax;
59 
60   f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
61   f[1] = (-Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s) + ctx->Pm) * ctx->omega_s / (2.0 * ctx->H);
62 
63   PetscCall(VecRestoreArrayRead(U, &u));
64   PetscCall(VecRestoreArray(F, &f));
65   PetscFunctionReturn(0);
66 }
67 
68 /*
69      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
70 */
71 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
72 {
73   PetscInt           rowcol[] = {0, 1};
74   PetscScalar        J[2][2], Pmax;
75   const PetscScalar *u;
76 
77   PetscFunctionBegin;
78   PetscCall(VecGetArrayRead(U, &u));
79   if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
80   else Pmax = ctx->Pmax;
81 
82   J[0][0] = 0;
83   J[0][1] = ctx->omega_b;
84   J[1][1] = -ctx->D * ctx->omega_s / (2.0 * ctx->H);
85   J[1][0] = -Pmax * PetscCosScalar(u[0]) * ctx->omega_s / (2.0 * ctx->H);
86 
87   PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
88   PetscCall(VecRestoreArrayRead(U, &u));
89 
90   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
91   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
92   if (A != B) {
93     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
94     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
95   }
96   PetscFunctionReturn(0);
97 }
98 
99 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx0)
100 {
101   PetscInt    row[] = {0, 1}, col[] = {0};
102   PetscScalar J[2][1];
103   AppCtx     *ctx = (AppCtx *)ctx0;
104 
105   PetscFunctionBeginUser;
106   J[0][0] = 0;
107   J[1][0] = ctx->omega_s / (2.0 * ctx->H);
108   PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
109   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
110   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
111   PetscFunctionReturn(0);
112 }
113 
114 static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
115 {
116   PetscScalar       *r;
117   const PetscScalar *u;
118 
119   PetscFunctionBegin;
120   PetscCall(VecGetArrayRead(U, &u));
121   PetscCall(VecGetArray(R, &r));
122   r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
123   PetscCall(VecRestoreArray(R, &r));
124   PetscCall(VecRestoreArrayRead(U, &u));
125   PetscFunctionReturn(0);
126 }
127 
128 static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
129 {
130   PetscScalar        ru[1];
131   const PetscScalar *u;
132   PetscInt           row[] = {0}, col[] = {0};
133 
134   PetscFunctionBegin;
135   PetscCall(VecGetArrayRead(U, &u));
136   ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
137   PetscCall(VecRestoreArrayRead(U, &u));
138   PetscCall(MatSetValues(DRDU, 1, row, 1, col, ru, INSERT_VALUES));
139   PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
140   PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
141   PetscFunctionReturn(0);
142 }
143 
144 static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, AppCtx *ctx)
145 {
146   PetscFunctionBegin;
147   PetscCall(MatZeroEntries(DRDP));
148   PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
149   PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
150   PetscFunctionReturn(0);
151 }
152 
153 PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
154 {
155   PetscScalar       *y, sensip;
156   const PetscScalar *x;
157 
158   PetscFunctionBegin;
159   PetscCall(VecGetArrayRead(lambda, &x));
160   PetscCall(VecGetArray(mu, &y));
161   sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
162   y[0]   = sensip;
163   PetscCall(VecRestoreArray(mu, &y));
164   PetscCall(VecRestoreArrayRead(lambda, &x));
165   PetscFunctionReturn(0);
166 }
167 
168 int main(int argc, char **argv)
169 {
170   Vec          p;
171   PetscScalar *x_ptr;
172   PetscMPIInt  size;
173   AppCtx       ctx;
174   Vec          lowerb, upperb;
175   Tao          tao;
176   KSP          ksp;
177   PC           pc;
178   Vec          U, lambda[1], mu[1];
179   Mat          A;    /* Jacobian matrix */
180   Mat          Jacp; /* Jacobian matrix */
181   Mat          DRDU, DRDP;
182   PetscInt     n = 2;
183   TS           quadts;
184 
185   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186      Initialize program
187      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188   PetscFunctionBeginUser;
189   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
190   PetscFunctionBeginUser;
191   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
192   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
193 
194   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195     Set runtime options
196     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
198   {
199     ctx.beta    = 2;
200     ctx.c       = PetscRealConstant(10000.0);
201     ctx.u_s     = PetscRealConstant(1.0);
202     ctx.omega_s = PetscRealConstant(1.0);
203     ctx.omega_b = PetscRealConstant(120.0) * PETSC_PI;
204     ctx.H       = PetscRealConstant(5.0);
205     PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
206     ctx.D = PetscRealConstant(5.0);
207     PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
208     ctx.E    = PetscRealConstant(1.1378);
209     ctx.V    = PetscRealConstant(1.0);
210     ctx.X    = PetscRealConstant(0.545);
211     ctx.Pmax = ctx.E * ctx.V / ctx.X;
212     PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
213     ctx.Pm = PetscRealConstant(1.0194);
214     PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
215     ctx.tf  = PetscRealConstant(0.1);
216     ctx.tcl = PetscRealConstant(0.2);
217     PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
218     PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
219   }
220   PetscOptionsEnd();
221 
222   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
223     Create necessary matrix and vectors
224     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
225   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
226   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
227   PetscCall(MatSetType(A, MATDENSE));
228   PetscCall(MatSetFromOptions(A));
229   PetscCall(MatSetUp(A));
230 
231   PetscCall(MatCreateVecs(A, &U, NULL));
232 
233   PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
234   PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
235   PetscCall(MatSetFromOptions(Jacp));
236   PetscCall(MatSetUp(Jacp));
237   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &DRDP));
238   PetscCall(MatSetUp(DRDP));
239   PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, &DRDU));
240   PetscCall(MatSetUp(DRDU));
241 
242   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243      Create timestepping solver context
244      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
245   PetscCall(TSCreate(PETSC_COMM_WORLD, &ctx.ts));
246   PetscCall(TSSetProblemType(ctx.ts, TS_NONLINEAR));
247   PetscCall(TSSetEquationType(ctx.ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
248   PetscCall(TSSetType(ctx.ts, TSRK));
249   PetscCall(TSSetRHSFunction(ctx.ts, NULL, (TSRHSFunction)RHSFunction, &ctx));
250   PetscCall(TSSetRHSJacobian(ctx.ts, A, A, (TSRHSJacobian)RHSJacobian, &ctx));
251   PetscCall(TSSetExactFinalTime(ctx.ts, TS_EXACTFINALTIME_MATCHSTEP));
252 
253   PetscCall(MatCreateVecs(A, &lambda[0], NULL));
254   PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
255   PetscCall(TSSetCostGradients(ctx.ts, 1, lambda, mu));
256   PetscCall(TSSetRHSJacobianP(ctx.ts, Jacp, RHSJacobianP, &ctx));
257 
258   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259      Set solver options
260    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
261   PetscCall(TSSetMaxTime(ctx.ts, PetscRealConstant(1.0)));
262   PetscCall(TSSetTimeStep(ctx.ts, PetscRealConstant(0.01)));
263   PetscCall(TSSetFromOptions(ctx.ts));
264 
265   PetscCall(TSGetTimeStep(ctx.ts, &ctx.dt)); /* save the stepsize */
266 
267   PetscCall(TSCreateQuadratureTS(ctx.ts, PETSC_TRUE, &quadts));
268   PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunction)CostIntegrand, &ctx));
269   PetscCall(TSSetRHSJacobian(quadts, DRDU, DRDU, (TSRHSJacobian)DRDUJacobianTranspose, &ctx));
270   PetscCall(TSSetRHSJacobianP(quadts, DRDP, (TSRHSJacobianP)DRDPJacobianTranspose, &ctx));
271   PetscCall(TSSetSolution(ctx.ts, U));
272 
273   /* Create TAO solver and set desired solution method */
274   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
275   PetscCall(TaoSetType(tao, TAOBLMVM));
276 
277   /*
278      Optimization starts
279   */
280   /* Set initial solution guess */
281   PetscCall(VecCreateSeq(PETSC_COMM_WORLD, 1, &p));
282   PetscCall(VecGetArray(p, &x_ptr));
283   x_ptr[0] = ctx.Pm;
284   PetscCall(VecRestoreArray(p, &x_ptr));
285 
286   PetscCall(TaoSetSolution(tao, p));
287   /* Set routine for function and gradient evaluation */
288   PetscCall(TaoSetObjective(tao, FormFunction, (void *)&ctx));
289   PetscCall(TaoSetGradient(tao, NULL, FormGradient, (void *)&ctx));
290 
291   /* Set bounds for the optimization */
292   PetscCall(VecDuplicate(p, &lowerb));
293   PetscCall(VecDuplicate(p, &upperb));
294   PetscCall(VecGetArray(lowerb, &x_ptr));
295   x_ptr[0] = 0.;
296   PetscCall(VecRestoreArray(lowerb, &x_ptr));
297   PetscCall(VecGetArray(upperb, &x_ptr));
298   x_ptr[0] = PetscRealConstant(1.1);
299   PetscCall(VecRestoreArray(upperb, &x_ptr));
300   PetscCall(TaoSetVariableBounds(tao, lowerb, upperb));
301 
302   /* Check for any TAO command line options */
303   PetscCall(TaoSetFromOptions(tao));
304   PetscCall(TaoGetKSP(tao, &ksp));
305   if (ksp) {
306     PetscCall(KSPGetPC(ksp, &pc));
307     PetscCall(PCSetType(pc, PCNONE));
308   }
309 
310   /* SOLVE THE APPLICATION */
311   PetscCall(TaoSolve(tao));
312 
313   PetscCall(VecView(p, PETSC_VIEWER_STDOUT_WORLD));
314   /* Free TAO data structures */
315   PetscCall(TaoDestroy(&tao));
316   PetscCall(VecDestroy(&p));
317   PetscCall(VecDestroy(&lowerb));
318   PetscCall(VecDestroy(&upperb));
319 
320   PetscCall(TSDestroy(&ctx.ts));
321   PetscCall(VecDestroy(&U));
322   PetscCall(MatDestroy(&A));
323   PetscCall(MatDestroy(&Jacp));
324   PetscCall(MatDestroy(&DRDU));
325   PetscCall(MatDestroy(&DRDP));
326   PetscCall(VecDestroy(&lambda[0]));
327   PetscCall(VecDestroy(&mu[0]));
328   PetscCall(PetscFinalize());
329   return 0;
330 }
331 
332 /* ------------------------------------------------------------------ */
333 /*
334    FormFunction - Evaluates the function
335 
336    Input Parameters:
337    tao - the Tao context
338    X   - the input vector
339    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
340 
341    Output Parameters:
342    f   - the newly evaluated function
343 */
344 PetscErrorCode FormFunction(Tao tao, Vec P, PetscReal *f, void *ctx0)
345 {
346   AppCtx      *ctx = (AppCtx *)ctx0;
347   TS           ts  = ctx->ts;
348   Vec          U; /* solution will be stored here */
349   PetscScalar *u;
350   PetscScalar *x_ptr;
351   Vec          q;
352 
353   PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
354   ctx->Pm = x_ptr[0];
355   PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
356 
357   /* reset time */
358   PetscCall(TSSetTime(ts, 0.0));
359   /* reset step counter, this is critical for adjoint solver */
360   PetscCall(TSSetStepNumber(ts, 0));
361   /* reset step size, the step size becomes negative after TSAdjointSolve */
362   PetscCall(TSSetTimeStep(ts, ctx->dt));
363   /* reinitialize the integral value */
364   PetscCall(TSGetCostIntegral(ts, &q));
365   PetscCall(VecSet(q, 0.0));
366 
367   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
368      Set initial conditions
369    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
370   PetscCall(TSGetSolution(ts, &U));
371   PetscCall(VecGetArray(U, &u));
372   u[0] = PetscAsinScalar(ctx->Pm / ctx->Pmax);
373   u[1] = PetscRealConstant(1.0);
374   PetscCall(VecRestoreArray(U, &u));
375 
376   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
377      Solve nonlinear system
378      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
379   PetscCall(TSSolve(ts, U));
380   PetscCall(TSGetCostIntegral(ts, &q));
381   PetscCall(VecGetArray(q, &x_ptr));
382   *f = -ctx->Pm + x_ptr[0];
383   PetscCall(VecRestoreArray(q, &x_ptr));
384   return 0;
385 }
386 
387 PetscErrorCode FormGradient(Tao tao, Vec P, Vec G, void *ctx0)
388 {
389   AppCtx      *ctx = (AppCtx *)ctx0;
390   TS           ts  = ctx->ts;
391   Vec          U; /* solution will be stored here */
392   PetscReal    ftime;
393   PetscInt     steps;
394   PetscScalar *u;
395   PetscScalar *x_ptr, *y_ptr;
396   Vec         *lambda, q, *mu;
397 
398   PetscCall(VecGetArrayRead(P, (const PetscScalar **)&x_ptr));
399   ctx->Pm = x_ptr[0];
400   PetscCall(VecRestoreArrayRead(P, (const PetscScalar **)&x_ptr));
401 
402   /* reset time */
403   PetscCall(TSSetTime(ts, 0.0));
404   /* reset step counter, this is critical for adjoint solver */
405   PetscCall(TSSetStepNumber(ts, 0));
406   /* reset step size, the step size becomes negative after TSAdjointSolve */
407   PetscCall(TSSetTimeStep(ts, ctx->dt));
408   /* reinitialize the integral value */
409   PetscCall(TSGetCostIntegral(ts, &q));
410   PetscCall(VecSet(q, 0.0));
411 
412   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413      Set initial conditions
414    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
415   PetscCall(TSGetSolution(ts, &U));
416   PetscCall(VecGetArray(U, &u));
417   u[0] = PetscAsinScalar(ctx->Pm / ctx->Pmax);
418   u[1] = PetscRealConstant(1.0);
419   PetscCall(VecRestoreArray(U, &u));
420 
421   /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */
422   PetscCall(TSSetSaveTrajectory(ts));
423   PetscCall(TSSetFromOptions(ts));
424 
425   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
426      Solve nonlinear system
427      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
428   PetscCall(TSSolve(ts, U));
429 
430   PetscCall(TSGetSolveTime(ts, &ftime));
431   PetscCall(TSGetStepNumber(ts, &steps));
432 
433   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434      Adjoint model starts here
435      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
436   PetscCall(TSGetCostGradients(ts, NULL, &lambda, &mu));
437   /*   Set initial conditions for the adjoint integration */
438   PetscCall(VecGetArray(lambda[0], &y_ptr));
439   y_ptr[0] = 0.0;
440   y_ptr[1] = 0.0;
441   PetscCall(VecRestoreArray(lambda[0], &y_ptr));
442   PetscCall(VecGetArray(mu[0], &x_ptr));
443   x_ptr[0] = PetscRealConstant(-1.0);
444   PetscCall(VecRestoreArray(mu[0], &x_ptr));
445 
446   PetscCall(TSAdjointSolve(ts));
447   PetscCall(TSGetCostIntegral(ts, &q));
448   PetscCall(ComputeSensiP(lambda[0], mu[0], ctx));
449   PetscCall(VecCopy(mu[0], G));
450   return 0;
451 }
452 
453 /*TEST
454 
455    build:
456       requires: !complex
457 
458    test:
459       args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason
460 
461    test:
462       suffix: 2
463       args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason -tao_test_gradient
464 
465 TEST*/
466