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; J[0][1] = ctx->omega_b; 83 J[1][1] = -ctx->D*ctx->omega_s/(2.0*ctx->H); J[1][0] = -Pmax*PetscCosScalar(u[0])*ctx->omega_s/(2.0*ctx->H); 84 85 PetscCall(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 86 PetscCall(VecRestoreArrayRead(U,&u)); 87 88 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 89 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 90 if (A != B) { 91 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 92 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 93 } 94 PetscFunctionReturn(0); 95 } 96 97 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0) 98 { 99 PetscInt row[] = {0,1},col[]={0}; 100 PetscScalar J[2][1]; 101 AppCtx *ctx=(AppCtx*)ctx0; 102 103 PetscFunctionBeginUser; 104 J[0][0] = 0; 105 J[1][0] = ctx->omega_s/(2.0*ctx->H); 106 PetscCall(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 107 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 108 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx) 113 { 114 PetscScalar *r; 115 const PetscScalar *u; 116 117 PetscFunctionBegin; 118 PetscCall(VecGetArrayRead(U,&u)); 119 PetscCall(VecGetArray(R,&r)); 120 r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta); 121 PetscCall(VecRestoreArray(R,&r)); 122 PetscCall(VecRestoreArrayRead(U,&u)); 123 PetscFunctionReturn(0); 124 } 125 126 static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx) 127 { 128 PetscScalar ru[1]; 129 const PetscScalar *u; 130 PetscInt row[] = {0},col[] = {0}; 131 132 PetscFunctionBegin; 133 PetscCall(VecGetArrayRead(U,&u)); 134 ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta-1); 135 PetscCall(VecRestoreArrayRead(U,&u)); 136 PetscCall(MatSetValues(DRDU,1,row,1,col,ru,INSERT_VALUES)); 137 PetscCall(MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY)); 138 PetscCall(MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY)); 139 PetscFunctionReturn(0); 140 } 141 142 static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx) 143 { 144 PetscFunctionBegin; 145 PetscCall(MatZeroEntries(DRDP)); 146 PetscCall(MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY)); 147 PetscCall(MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY)); 148 PetscFunctionReturn(0); 149 } 150 151 PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx) 152 { 153 PetscScalar *y,sensip; 154 const PetscScalar *x; 155 156 PetscFunctionBegin; 157 PetscCall(VecGetArrayRead(lambda,&x)); 158 PetscCall(VecGetArray(mu,&y)); 159 sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0]; 160 y[0] = sensip; 161 PetscCall(VecRestoreArray(mu,&y)); 162 PetscCall(VecRestoreArrayRead(lambda,&x)); 163 PetscFunctionReturn(0); 164 } 165 166 int main(int argc,char **argv) 167 { 168 Vec p; 169 PetscScalar *x_ptr; 170 PetscMPIInt size; 171 AppCtx ctx; 172 Vec lowerb,upperb; 173 Tao tao; 174 KSP ksp; 175 PC pc; 176 Vec U,lambda[1],mu[1]; 177 Mat A; /* Jacobian matrix */ 178 Mat Jacp; /* Jacobian matrix */ 179 Mat DRDU,DRDP; 180 PetscInt n = 2; 181 TS quadts; 182 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Initialize program 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 187 PetscFunctionBeginUser; 188 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 189 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 190 191 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 192 Set runtime options 193 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 194 PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options",""); 195 { 196 ctx.beta = 2; 197 ctx.c = PetscRealConstant(10000.0); 198 ctx.u_s = PetscRealConstant(1.0); 199 ctx.omega_s = PetscRealConstant(1.0); 200 ctx.omega_b = PetscRealConstant(120.0)*PETSC_PI; 201 ctx.H = PetscRealConstant(5.0); 202 PetscCall(PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL)); 203 ctx.D = PetscRealConstant(5.0); 204 PetscCall(PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL)); 205 ctx.E = PetscRealConstant(1.1378); 206 ctx.V = PetscRealConstant(1.0); 207 ctx.X = PetscRealConstant(0.545); 208 ctx.Pmax = ctx.E*ctx.V/ctx.X; 209 PetscCall(PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL)); 210 ctx.Pm = PetscRealConstant(1.0194); 211 PetscCall(PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL)); 212 ctx.tf = PetscRealConstant(0.1); 213 ctx.tcl = PetscRealConstant(0.2); 214 PetscCall(PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL)); 215 PetscCall(PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL)); 216 217 } 218 PetscOptionsEnd(); 219 220 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 221 Create necessary matrix and vectors 222 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 223 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 224 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 225 PetscCall(MatSetType(A,MATDENSE)); 226 PetscCall(MatSetFromOptions(A)); 227 PetscCall(MatSetUp(A)); 228 229 PetscCall(MatCreateVecs(A,&U,NULL)); 230 231 PetscCall(MatCreate(PETSC_COMM_WORLD,&Jacp)); 232 PetscCall(MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1)); 233 PetscCall(MatSetFromOptions(Jacp)); 234 PetscCall(MatSetUp(Jacp)); 235 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&DRDP)); 236 PetscCall(MatSetUp(DRDP)); 237 PetscCall(MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,2,NULL,&DRDU)); 238 PetscCall(MatSetUp(DRDU)); 239 240 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241 Create timestepping solver context 242 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 243 PetscCall(TSCreate(PETSC_COMM_WORLD,&ctx.ts)); 244 PetscCall(TSSetProblemType(ctx.ts,TS_NONLINEAR)); 245 PetscCall(TSSetEquationType(ctx.ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 246 PetscCall(TSSetType(ctx.ts,TSRK)); 247 PetscCall(TSSetRHSFunction(ctx.ts,NULL,(TSRHSFunction)RHSFunction,&ctx)); 248 PetscCall(TSSetRHSJacobian(ctx.ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx)); 249 PetscCall(TSSetExactFinalTime(ctx.ts,TS_EXACTFINALTIME_MATCHSTEP)); 250 251 PetscCall(MatCreateVecs(A,&lambda[0],NULL)); 252 PetscCall(MatCreateVecs(Jacp,&mu[0],NULL)); 253 PetscCall(TSSetCostGradients(ctx.ts,1,lambda,mu)); 254 PetscCall(TSSetRHSJacobianP(ctx.ts,Jacp,RHSJacobianP,&ctx)); 255 256 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 257 Set solver options 258 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 259 PetscCall(TSSetMaxTime(ctx.ts,PetscRealConstant(1.0))); 260 PetscCall(TSSetTimeStep(ctx.ts,PetscRealConstant(0.01))); 261 PetscCall(TSSetFromOptions(ctx.ts)); 262 263 PetscCall(TSGetTimeStep(ctx.ts,&ctx.dt)); /* save the stepsize */ 264 265 PetscCall(TSCreateQuadratureTS(ctx.ts,PETSC_TRUE,&quadts)); 266 PetscCall(TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx)); 267 PetscCall(TSSetRHSJacobian(quadts,DRDU,DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx)); 268 PetscCall(TSSetRHSJacobianP(quadts,DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&ctx)); 269 PetscCall(TSSetSolution(ctx.ts,U)); 270 271 /* Create TAO solver and set desired solution method */ 272 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 273 PetscCall(TaoSetType(tao,TAOBLMVM)); 274 275 /* 276 Optimization starts 277 */ 278 /* Set initial solution guess */ 279 PetscCall(VecCreateSeq(PETSC_COMM_WORLD,1,&p)); 280 PetscCall(VecGetArray(p,&x_ptr)); 281 x_ptr[0] = ctx.Pm; 282 PetscCall(VecRestoreArray(p,&x_ptr)); 283 284 PetscCall(TaoSetSolution(tao,p)); 285 /* Set routine for function and gradient evaluation */ 286 PetscCall(TaoSetObjective(tao,FormFunction,(void *)&ctx)); 287 PetscCall(TaoSetGradient(tao,NULL,FormGradient,(void *)&ctx)); 288 289 /* Set bounds for the optimization */ 290 PetscCall(VecDuplicate(p,&lowerb)); 291 PetscCall(VecDuplicate(p,&upperb)); 292 PetscCall(VecGetArray(lowerb,&x_ptr)); 293 x_ptr[0] = 0.; 294 PetscCall(VecRestoreArray(lowerb,&x_ptr)); 295 PetscCall(VecGetArray(upperb,&x_ptr)); 296 x_ptr[0] = PetscRealConstant(1.1); 297 PetscCall(VecRestoreArray(upperb,&x_ptr)); 298 PetscCall(TaoSetVariableBounds(tao,lowerb,upperb)); 299 300 /* Check for any TAO command line options */ 301 PetscCall(TaoSetFromOptions(tao)); 302 PetscCall(TaoGetKSP(tao,&ksp)); 303 if (ksp) { 304 PetscCall(KSPGetPC(ksp,&pc)); 305 PetscCall(PCSetType(pc,PCNONE)); 306 } 307 308 /* SOLVE THE APPLICATION */ 309 PetscCall(TaoSolve(tao)); 310 311 PetscCall(VecView(p,PETSC_VIEWER_STDOUT_WORLD)); 312 /* Free TAO data structures */ 313 PetscCall(TaoDestroy(&tao)); 314 PetscCall(VecDestroy(&p)); 315 PetscCall(VecDestroy(&lowerb)); 316 PetscCall(VecDestroy(&upperb)); 317 318 PetscCall(TSDestroy(&ctx.ts)); 319 PetscCall(VecDestroy(&U)); 320 PetscCall(MatDestroy(&A)); 321 PetscCall(MatDestroy(&Jacp)); 322 PetscCall(MatDestroy(&DRDU)); 323 PetscCall(MatDestroy(&DRDP)); 324 PetscCall(VecDestroy(&lambda[0])); 325 PetscCall(VecDestroy(&mu[0])); 326 PetscCall(PetscFinalize()); 327 return 0; 328 } 329 330 /* ------------------------------------------------------------------ */ 331 /* 332 FormFunction - Evaluates the function 333 334 Input Parameters: 335 tao - the Tao context 336 X - the input vector 337 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 338 339 Output Parameters: 340 f - the newly evaluated function 341 */ 342 PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0) 343 { 344 AppCtx *ctx = (AppCtx*)ctx0; 345 TS ts = ctx->ts; 346 Vec U; /* solution will be stored here */ 347 PetscScalar *u; 348 PetscScalar *x_ptr; 349 Vec q; 350 351 PetscCall(VecGetArrayRead(P,(const PetscScalar**)&x_ptr)); 352 ctx->Pm = x_ptr[0]; 353 PetscCall(VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr)); 354 355 /* reset time */ 356 PetscCall(TSSetTime(ts,0.0)); 357 /* reset step counter, this is critical for adjoint solver */ 358 PetscCall(TSSetStepNumber(ts,0)); 359 /* reset step size, the step size becomes negative after TSAdjointSolve */ 360 PetscCall(TSSetTimeStep(ts,ctx->dt)); 361 /* reinitialize the integral value */ 362 PetscCall(TSGetCostIntegral(ts,&q)); 363 PetscCall(VecSet(q,0.0)); 364 365 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 366 Set initial conditions 367 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 368 PetscCall(TSGetSolution(ts,&U)); 369 PetscCall(VecGetArray(U,&u)); 370 u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax); 371 u[1] = PetscRealConstant(1.0); 372 PetscCall(VecRestoreArray(U,&u)); 373 374 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 375 Solve nonlinear system 376 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 377 PetscCall(TSSolve(ts,U)); 378 PetscCall(TSGetCostIntegral(ts,&q)); 379 PetscCall(VecGetArray(q,&x_ptr)); 380 *f = -ctx->Pm + x_ptr[0]; 381 PetscCall(VecRestoreArray(q,&x_ptr)); 382 return 0; 383 } 384 385 PetscErrorCode FormGradient(Tao tao,Vec P,Vec G,void *ctx0) 386 { 387 AppCtx *ctx = (AppCtx*)ctx0; 388 TS ts = ctx->ts; 389 Vec U; /* solution will be stored here */ 390 PetscReal ftime; 391 PetscInt steps; 392 PetscScalar *u; 393 PetscScalar *x_ptr,*y_ptr; 394 Vec *lambda,q,*mu; 395 396 PetscCall(VecGetArrayRead(P,(const PetscScalar**)&x_ptr)); 397 ctx->Pm = x_ptr[0]; 398 PetscCall(VecRestoreArrayRead(P,(const PetscScalar**)&x_ptr)); 399 400 /* reset time */ 401 PetscCall(TSSetTime(ts,0.0)); 402 /* reset step counter, this is critical for adjoint solver */ 403 PetscCall(TSSetStepNumber(ts,0)); 404 /* reset step size, the step size becomes negative after TSAdjointSolve */ 405 PetscCall(TSSetTimeStep(ts,ctx->dt)); 406 /* reinitialize the integral value */ 407 PetscCall(TSGetCostIntegral(ts,&q)); 408 PetscCall(VecSet(q,0.0)); 409 410 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 411 Set initial conditions 412 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 413 PetscCall(TSGetSolution(ts,&U)); 414 PetscCall(VecGetArray(U,&u)); 415 u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax); 416 u[1] = PetscRealConstant(1.0); 417 PetscCall(VecRestoreArray(U,&u)); 418 419 /* Set up to save trajectory before TSSetFromOptions() so that TSTrajectory options can be captured */ 420 PetscCall(TSSetSaveTrajectory(ts)); 421 PetscCall(TSSetFromOptions(ts)); 422 423 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 424 Solve nonlinear system 425 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 426 PetscCall(TSSolve(ts,U)); 427 428 PetscCall(TSGetSolveTime(ts,&ftime)); 429 PetscCall(TSGetStepNumber(ts,&steps)); 430 431 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 432 Adjoint model starts here 433 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 434 PetscCall(TSGetCostGradients(ts,NULL,&lambda,&mu)); 435 /* Set initial conditions for the adjoint integration */ 436 PetscCall(VecGetArray(lambda[0],&y_ptr)); 437 y_ptr[0] = 0.0; y_ptr[1] = 0.0; 438 PetscCall(VecRestoreArray(lambda[0],&y_ptr)); 439 PetscCall(VecGetArray(mu[0],&x_ptr)); 440 x_ptr[0] = PetscRealConstant(-1.0); 441 PetscCall(VecRestoreArray(mu[0],&x_ptr)); 442 443 PetscCall(TSAdjointSolve(ts)); 444 PetscCall(TSGetCostIntegral(ts,&q)); 445 PetscCall(ComputeSensiP(lambda[0],mu[0],ctx)); 446 PetscCall(VecCopy(mu[0],G)); 447 return 0; 448 } 449 450 /*TEST 451 452 build: 453 requires: !complex 454 455 test: 456 args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason 457 458 test: 459 suffix: 2 460 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 461 462 TEST*/ 463