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