Actual source code: ex9opt.c
petsc-3.11.3 2019-06-26
2: static char help[] = "Basic equation for generator stability analysis.\n" ;
\begin{eqnarray}
\frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
\frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
\end{eqnarray}
Ensemble of initial conditions
./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
Fault at .1 seconds
./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
Initial conditions same as when fault is ended
./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
25: /*
26: Include "petscts.h" so that we can use TS solvers. Note that this
27: file automatically includes:
28: petscsys.h - base PETSc routines petscvec.h - vectors
29: petscmat.h - matrices
30: petscis.h - index sets petscksp.h - Krylov subspace methods
31: petscviewer.h - viewers petscpc.h - preconditioners
32: petscksp.h - linear solvers
33: */
35: #include <petsctao.h>
36: #include <petscts.h>
38: typedef struct {
39: PetscScalar H,D,omega_b,omega_s,Pmax,Pm,E,V,X,u_s,c;
40: PetscInt beta;
41: PetscReal tf,tcl;
42: } AppCtx;
44: PetscErrorCode FormFunction(Tao ,Vec ,PetscReal *,void*) ;
45: PetscErrorCode FormGradient(Tao ,Vec ,Vec ,void*) ;
47: /*
48: Defines the ODE passed to the ODE solver
49: */
50: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
51: {
52: PetscErrorCode ierr;
53: PetscScalar *f,Pmax;
54: const PetscScalar *u;
57: /* The next three lines allow us to access the entries of the vectors directly */
58: VecGetArrayRead (U,&u);
59: VecGetArray (F,&f);
60: 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 */
61: else Pmax = ctx->Pmax;
63: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
64: f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H);
66: VecRestoreArrayRead (U,&u);
67: VecRestoreArray (F,&f);
68: return (0);
69: }
71: /*
72: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian () for the meaning of a and the Jacobian.
73: */
74: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
75: {
76: PetscErrorCode ierr;
77: PetscInt rowcol[] = {0,1};
78: PetscScalar J[2][2],Pmax;
79: const PetscScalar *u;
82: VecGetArrayRead (U,&u);
83: 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 */
84: else Pmax = ctx->Pmax;
86: J[0][0] = 0; J[0][1] = ctx->omega_b;
87: 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);
89: MatSetValues (B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES );
90: VecRestoreArrayRead (U,&u);
92: MatAssemblyBegin (A,MAT_FINAL_ASSEMBLY );
93: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY );
94: if (A != B) {
95: MatAssemblyBegin (B,MAT_FINAL_ASSEMBLY );
96: MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY );
97: }
98: return (0);
99: }
101: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
102: {
104: PetscInt row[] = {0,1},col[]={0};
105: PetscScalar J[2][1];
106: AppCtx *ctx=(AppCtx*)ctx0;
109: J[0][0] = 0;
110: J[1][0] = ctx->omega_s/(2.0*ctx->H);
111: MatSetValues (A,2,row,1,col,&J[0][0],INSERT_VALUES );
112: MatAssemblyBegin (A,MAT_FINAL_ASSEMBLY );
113: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY );
114: return (0);
115: }
117: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
118: {
119: PetscErrorCode ierr;
120: PetscScalar *r;
121: const PetscScalar *u;
124: VecGetArrayRead (U,&u);
125: VecGetArray (R,&r);
126: r[0] = ctx->c*PetscPowScalarInt(PetscMax (0., u[0]-ctx->u_s),ctx->beta);
127: VecRestoreArray (R,&r);
128: VecRestoreArrayRead (U,&u);
129: return (0);
130: }
132: static PetscErrorCode DRDYFunction(TS ts,PetscReal t,Vec U,Vec *drdy,AppCtx *ctx)
133: {
134: PetscErrorCode ierr;
135: PetscScalar *ry;
136: const PetscScalar *u;
139: VecGetArrayRead (U,&u);
140: VecGetArray (drdy[0],&ry);
141: ry[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax (0., u[0]-ctx->u_s),ctx->beta-1);
142: VecRestoreArray (drdy[0],&ry);
143: VecRestoreArrayRead (U,&u);
144: return (0);
145: }
147: static PetscErrorCode DRDPFunction(TS ts,PetscReal t,Vec U,Vec *drdp,AppCtx *ctx)
148: {
150: PetscScalar *rp;
153: VecGetArray (drdp[0],&rp);
154: rp[0] = 0.;
155: VecRestoreArray (drdp[0],&rp);
156: return (0);
157: }
159: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
160: {
161: PetscErrorCode ierr;
162: PetscScalar *y,sensip;
163: const PetscScalar *x;
166: VecGetArrayRead (lambda,&x);
167: VecGetArray (mu,&y);
168: sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
169: /*PetscPrintf (PETSC_COMM_WORLD ,"\n sensitivity wrt parameter pm: %g \n",(double)sensip);*/
170: y[0] = sensip;
171: VecRestoreArray (mu,&y);
172: VecRestoreArrayRead (lambda,&x);
173: return (0);
174: }
176: int main(int argc,char **argv)
177: {
178: Vec p;
179: PetscScalar *x_ptr;
180: PetscErrorCode ierr;
181: PetscMPIInt size;
182: AppCtx ctx;
183: Vec lowerb,upperb;
184: Tao tao;
185: KSP ksp;
186: PC pc;
188: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
189: Initialize program
190: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
191: PetscInitialize (&argc,&argv,NULL,help);if (ierr) return ierr;
193: MPI_Comm_size (PETSC_COMM_WORLD ,&size);
194: if (size != 1) SETERRQ (PETSC_COMM_SELF ,1,"This is a uniprocessor example only!" );
196: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: Set runtime options
198: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199: PetscOptionsBegin (PETSC_COMM_WORLD ,NULL,"Swing equation options" ,"" );
200: {
201: ctx.beta = 2;
202: ctx.c = 10000.0;
203: ctx.u_s = 1.0;
204: ctx.omega_s = 1.0;
205: ctx.omega_b = 120.0*PETSC_PI;
206: ctx.H = 5.0;
207: PetscOptionsScalar ("-Inertia" ,"" ,"" ,ctx.H,&ctx.H,NULL);
208: ctx.D = 5.0;
209: PetscOptionsScalar ("-D" ,"" ,"" ,ctx.D,&ctx.D,NULL);
210: #if defined(PETSC_USE_REAL___FLOAT128)
211: ctx.E = 1.1378q;
212: #else
213: ctx.E = 1.1378;
214: #endif
215: ctx.V = 1.0;
216: #if defined(PETSC_USE_REAL___FLOAT128)
217: ctx.X = 0.545q;
218: #else
219: ctx.X = 0.545;
220: #endif
221: ctx.Pmax = ctx.E*ctx.V/ctx.X;;
222: PetscOptionsScalar ("-Pmax" ,"" ,"" ,ctx.Pmax,&ctx.Pmax,NULL);
223: #if defined(PETSC_USE_REAL___FLOAT128)
224: ctx.Pm = 1.0194q;
225: #else
226: ctx.Pm = 1.0194;
227: #endif
228: PetscOptionsScalar ("-Pm" ,"" ,"" ,ctx.Pm,&ctx.Pm,NULL);
229: #if defined(PETSC_USE_REAL___FLOAT128)
230: ctx.tf = 0.1q;
231: ctx.tcl = 0.2q;
232: #else
233: ctx.tf = 0.1;
234: ctx.tcl = 0.2;
235: #endif
236: PetscOptionsReal ("-tf" ,"Time to start fault" ,"" ,ctx.tf,&ctx.tf,NULL);
237: PetscOptionsReal ("-tcl" ,"Time to end fault" ,"" ,ctx.tcl,&ctx.tcl,NULL);
239: }
240: PetscOptionsEnd ();
242: /* Create TAO solver and set desired solution method */
243: TaoCreate (PETSC_COMM_WORLD ,&tao);
244: TaoSetType (tao,TAOBLMVM );
246: /*
247: Optimization starts
248: */
249: /* Set initial solution guess */
250: VecCreateSeq (PETSC_COMM_WORLD ,1,&p);
251: VecGetArray (p,&x_ptr);
252: x_ptr[0] = ctx.Pm;
253: VecRestoreArray (p,&x_ptr);
255: TaoSetInitialVector (tao,p);
256: /* Set routine for function and gradient evaluation */
257: TaoSetObjectiveRoutine (tao,FormFunction,(void *)&ctx);
258: TaoSetGradientRoutine (tao,FormGradient,(void *)&ctx);
260: /* Set bounds for the optimization */
261: VecDuplicate (p,&lowerb);
262: VecDuplicate (p,&upperb);
263: VecGetArray (lowerb,&x_ptr);
264: x_ptr[0] = 0.;
265: VecRestoreArray (lowerb,&x_ptr);
266: VecGetArray (upperb,&x_ptr);
267: #if defined(PETSC_USE_REAL___FLOAT128)
268: x_ptr[0] = 1.1q;
269: #else
270: x_ptr[0] = 1.1;
271: #endif
272: VecRestoreArray (upperb,&x_ptr);
273: TaoSetVariableBounds (tao,lowerb,upperb);
275: /* Check for any TAO command line options */
276: TaoSetFromOptions (tao);
277: TaoGetKSP (tao,&ksp);
278: if (ksp) {
279: KSPGetPC (ksp,&pc);
280: PCSetType (pc,PCNONE );
281: }
283: /* SOLVE THE APPLICATION */
284: TaoSolve (tao);
286: VecView (p,PETSC_VIEWER_STDOUT_WORLD );
287: /* Free TAO data structures */
288: TaoDestroy (&tao);
289: VecDestroy (&p);
290: VecDestroy (&lowerb);
291: VecDestroy (&upperb);
292: PetscFinalize ();
293: return ierr;
294: }
296: /* ------------------------------------------------------------------ */
297: /*
298: FormFunction - Evaluates the function
300: Input Parameters:
301: tao - the Tao context
302: X - the input vector
303: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine ()
305: Output Parameters:
306: f - the newly evaluated function
307: */
308: PetscErrorCode FormFunction(Tao tao,Vec P,PetscReal *f,void *ctx0)
309: {
310: AppCtx *ctx = (AppCtx*)ctx0;
311: TS ts;
313: Vec U; /* solution will be stored here */
314: Mat A; /* Jacobian matrix */
316: PetscInt n = 2;
317: PetscScalar *u;
318: PetscScalar *x_ptr;
319: Vec lambda[1],q,mu[1];
321: VecGetArrayRead (P,(const PetscScalar **)&x_ptr);
322: ctx->Pm = x_ptr[0];
323: VecRestoreArrayRead (P,(const PetscScalar **)&x_ptr);
324: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
325: Create necessary matrix and vectors
326: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
327: MatCreate (PETSC_COMM_WORLD ,&A);
328: MatSetSizes (A,n,n,PETSC_DETERMINE ,PETSC_DETERMINE );
329: MatSetType (A,MATDENSE );
330: MatSetFromOptions (A);
331: MatSetUp (A);
333: MatCreateVecs (A,&U,NULL);
334: MatCreateVecs (A,&lambda[0],NULL);
335: MatCreateVecs (A,&mu[0],NULL);
337: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
338: Create timestepping solver context
339: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
340: TSCreate (PETSC_COMM_WORLD ,&ts);
341: TSSetProblemType (ts,TS_NONLINEAR );
342: TSSetType (ts,TSRK );
343: TSSetRHSFunction (ts,NULL,(TSRHSFunction)RHSFunction,ctx);
344: TSSetRHSJacobian (ts,A,A,(TSRHSJacobian)RHSJacobian,ctx);
345: TSSetExactFinalTime (ts,TS_EXACTFINALTIME_MATCHSTEP );
347: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: Set initial conditions
349: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350: VecGetArray (U,&u);
351: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
352: u[1] = 1.0;
353: VecRestoreArray (U,&u);
354: TSSetSolution (ts,U);
356: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
357: Save trajectory of solution so that TSAdjointSolve () may be used
358: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
359: TSSetSaveTrajectory (ts);
361: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
362: Set solver options
363: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
364: TSSetMaxTime (ts,1.0);
365: #if defined(PETSC_USE_REAL___FLOAT128)
366: TSSetTimeStep (ts,.01q);
367: #else
368: TSSetTimeStep (ts,.01);
369: #endif
370: TSSetFromOptions (ts);
372: TSSetCostGradients (ts,1,lambda,mu);
373: TSSetCostIntegrand (ts,1,NULL,(PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec ,void*))CostIntegrand,
374: (PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec *,void*))DRDYFunction,
375: (PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec *,void*))DRDPFunction,PETSC_TRUE ,ctx);
377: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
378: Solve nonlinear system
379: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
380: TSSolve (ts,U);
381: TSGetCostIntegral (ts,&q);
382: VecGetArray (q,&x_ptr);
383: *f = -ctx->Pm + x_ptr[0];
384: VecRestoreArray (q,&x_ptr);
386: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
387: Free work space. All PETSc objects should be destroyed when they are no longer needed.
388: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
389: MatDestroy (&A);
390: VecDestroy (&U);
391: VecDestroy (&lambda[0]);
392: VecDestroy (&mu[0]);
393: TSDestroy (&ts);
395: return 0;
396: }
399: PetscErrorCode FormGradient(Tao tao,Vec P,Vec G,void *ctx0)
400: {
401: AppCtx *ctx = (AppCtx*)ctx0;
402: TS ts;
404: Vec U; /* solution will be stored here */
405: Mat A; /* Jacobian matrix */
406: Mat Jacp; /* Jacobian matrix */
408: PetscInt n = 2;
409: PetscReal ftime;
410: PetscInt steps;
411: PetscScalar *u;
412: PetscScalar *x_ptr,*y_ptr;
413: Vec lambda[1],q,mu[1];
415: VecGetArrayRead (P,(const PetscScalar **)&x_ptr);
416: ctx->Pm = x_ptr[0];
417: VecRestoreArrayRead (P,(const PetscScalar **)&x_ptr);
419: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
420: Create necessary matrix and vectors
421: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
422: MatCreate (PETSC_COMM_WORLD ,&A);
423: MatSetSizes (A,n,n,PETSC_DETERMINE ,PETSC_DETERMINE );
424: MatSetType (A,MATDENSE );
425: MatSetFromOptions (A);
426: MatSetUp (A);
428: MatCreateVecs (A,&U,NULL);
430: MatCreate (PETSC_COMM_WORLD ,&Jacp);
431: MatSetSizes (Jacp,PETSC_DECIDE ,PETSC_DECIDE ,2,1);
432: MatSetFromOptions (Jacp);
433: MatSetUp (Jacp);
435: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436: Create timestepping solver context
437: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
438: TSCreate (PETSC_COMM_WORLD ,&ts);
439: TSSetProblemType (ts,TS_NONLINEAR );
440: TSSetType (ts,TSRK );
441: TSSetRHSFunction (ts,NULL,(TSRHSFunction)RHSFunction,ctx);
442: TSSetRHSJacobian (ts,A,A,(TSRHSJacobian)RHSJacobian,ctx);
443: TSSetExactFinalTime (ts,TS_EXACTFINALTIME_MATCHSTEP );
445: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
446: Set initial conditions
447: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
448: VecGetArray (U,&u);
449: u[0] = PetscAsinScalar(ctx->Pm/ctx->Pmax);
450: u[1] = 1.0;
451: VecRestoreArray (U,&u);
452: TSSetSolution (ts,U);
454: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
455: Save trajectory of solution so that TSAdjointSolve () may be used
456: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
457: TSSetSaveTrajectory (ts);
459: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
460: Set solver options
461: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
462: TSSetMaxTime (ts,1.0);
463: #if defined(PETSC_USE_REAL___FLOAT128)
464: TSSetTimeStep (ts,.01q);
465: #else
466: TSSetTimeStep (ts,.01);
467: #endif
468: TSSetFromOptions (ts);
470: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
471: Solve nonlinear system
472: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
473: TSSolve (ts,U);
475: TSGetSolveTime (ts,&ftime);
476: TSGetStepNumber (ts,&steps);
478: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
479: Adjoint model starts here
480: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
481: MatCreateVecs (A,&lambda[0],NULL);
482: /* Set initial conditions for the adjoint integration */
483: VecGetArray (lambda[0],&y_ptr);
484: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
485: VecRestoreArray (lambda[0],&y_ptr);
487: MatCreateVecs (Jacp,&mu[0],NULL);
488: VecGetArray (mu[0],&x_ptr);
489: x_ptr[0] = -1.0;
490: VecRestoreArray (mu[0],&x_ptr);
491: TSSetCostGradients (ts,1,lambda,mu);
493: TSSetRHSJacobianP (ts,Jacp,RHSJacobianP,ctx);
495: TSSetCostIntegrand (ts,1,NULL,(PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec ,void*))CostIntegrand,
496: (PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec *,void*))DRDYFunction,
497: (PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec *,void*))DRDPFunction,PETSC_FALSE ,ctx);
499: TSAdjointSolve (ts);
500: TSGetCostIntegral (ts,&q);
501: ComputeSensiP(lambda[0],mu[0],ctx);
503: VecCopy (mu[0],G);
505: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
506: Free work space. All PETSc objects should be destroyed when they are no longer needed.
507: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
508: MatDestroy (&A);
509: MatDestroy (&Jacp);
510: VecDestroy (&U);
511: VecDestroy (&lambda[0]);
512: VecDestroy (&mu[0]);
513: TSDestroy (&ts);
515: return 0;
516: }
519: /*TEST
521: build:
522: requires: !complex
524: test:
525: args: -viewer_binary_skip_info -ts_adapt_type none -tao_monitor -tao_gatol 0.0 -tao_grtol 1.e-3 -tao_converged_reason
527: test:
528: suffix: 2
529: 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
531: TEST*/