Actual source code: ex9adj.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
./ex9 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
Fault at .1 seconds
./ex9 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
Initial conditions same as when fault is ended
./ex9 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -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 <petscts.h>
37: typedef struct {
38: PetscScalar H,D,omega_b,omega_s,Pmax,Pm,E,V,X,u_s,c;
39: PetscInt beta;
40: PetscReal tf,tcl;
41: } AppCtx;
43: PetscErrorCode PostStepFunction(TS ts)
44: {
45: PetscErrorCode ierr;
46: Vec U;
47: PetscReal t;
48: const PetscScalar *u;
51: TSGetTime (ts,&t);
52: TSGetSolution (ts,&U);
53: VecGetArrayRead (U,&u);
54: PetscPrintf (PETSC_COMM_SELF ,"delta(%3.2f) = %8.7f\n" ,(double)t,(double)u[0]);
55: VecRestoreArrayRead (U,&u);
56:
57: return (0);
58: }
60: /*
61: Defines the ODE passed to the ODE solver
62: */
63: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
64: {
65: PetscErrorCode ierr;
66: PetscScalar *f,Pmax;
67: const PetscScalar *u;
70: /* The next three lines allow us to access the entries of the vectors directly */
71: VecGetArrayRead (U,&u);
72: VecGetArray (F,&f);
73: 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 */
74: else Pmax = ctx->Pmax;
75:
76: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
77: f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H);
79: VecRestoreArrayRead (U,&u);
80: VecRestoreArray (F,&f);
81: return (0);
82: }
84: /*
85: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian () for the meaning of a and the Jacobian.
86: */
87: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
88: {
89: PetscErrorCode ierr;
90: PetscInt rowcol[] = {0,1};
91: PetscScalar J[2][2],Pmax;
92: const PetscScalar *u;
95: VecGetArrayRead (U,&u);
96: 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 */
97: else Pmax = ctx->Pmax;
99: J[0][0] = 0; J[0][1] = ctx->omega_b;
100: 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);
102: MatSetValues (A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES );
103: VecRestoreArrayRead (U,&u);
105: MatAssemblyBegin (A,MAT_FINAL_ASSEMBLY );
106: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY );
107: if (A != B) {
108: MatAssemblyBegin (B,MAT_FINAL_ASSEMBLY );
109: MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY );
110: }
111: return (0);
112: }
114: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
115: {
117: PetscInt row[] = {0,1},col[]={0};
118: PetscScalar J[2][1];
119: AppCtx *ctx=(AppCtx*)ctx0;
122: J[0][0] = 0;
123: J[1][0] = ctx->omega_s/(2.0*ctx->H);
124: MatSetValues (A,2,row,1,col,&J[0][0],INSERT_VALUES );
125: MatAssemblyBegin (A,MAT_FINAL_ASSEMBLY );
126: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY );
127: return (0);
128: }
130: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
131: {
132: PetscErrorCode ierr;
133: PetscScalar *r;
134: const PetscScalar *u;
137: VecGetArrayRead (U,&u);
138: VecGetArray (R,&r);
139: r[0] = ctx->c*PetscPowScalarInt(PetscMax (0., u[0]-ctx->u_s),ctx->beta);
140: VecRestoreArray (R,&r);
141: VecRestoreArrayRead (U,&u);
142: return (0);
143: }
145: static PetscErrorCode DRDYFunction(TS ts,PetscReal t,Vec U,Vec *drdy,AppCtx *ctx)
146: {
147: PetscErrorCode ierr;
148: PetscScalar *ry;
149: const PetscScalar *u;
152: VecGetArrayRead (U,&u);
153: VecGetArray (drdy[0],&ry);
154: ry[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax (0., u[0]-ctx->u_s),ctx->beta-1);
155: VecRestoreArray (drdy[0],&ry);
156: VecRestoreArrayRead (U,&u);
157: return (0);
158: }
160: static PetscErrorCode DRDPFunction(TS ts,PetscReal t,Vec U,Vec *drdp,AppCtx *ctx)
161: {
162: PetscErrorCode ierr;
163: PetscScalar *rp;
164: const PetscScalar *u;
167: VecGetArrayRead (U,&u);
168: VecGetArray (drdp[0],&rp);
169: rp[0] = 0.;
170: VecRestoreArray (drdp[0],&rp);
171: VecGetArrayRead (U,&u);
172: return (0);
173: }
175: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
176: {
177: PetscErrorCode ierr;
178: PetscScalar sensip;
179: const PetscScalar *x,*y;
180:
182: VecGetArrayRead (lambda,&x);
183: VecGetArrayRead (mu,&y);
184: sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
185: PetscPrintf (PETSC_COMM_WORLD ,"\n sensitivity wrt parameter pm: %.7f \n" ,(double)sensip);
186: VecRestoreArrayRead (lambda,&x);
187: VecRestoreArrayRead (mu,&y);
188: return (0);
189: }
191: int main(int argc,char **argv)
192: {
193: TS ts; /* ODE integrator */
194: Vec U; /* solution will be stored here */
195: Mat A; /* Jacobian matrix */
196: Mat Jacp; /* Jacobian matrix */
198: PetscMPIInt size;
199: PetscInt n = 2;
200: AppCtx ctx;
201: PetscScalar *u;
202: PetscReal du[2] = {0.0,0.0};
203: PetscBool ensemble = PETSC_FALSE ,flg1,flg2;
204: PetscReal ftime;
205: PetscInt steps;
206: PetscScalar *x_ptr,*y_ptr;
207: Vec lambda[1],q,mu[1];
209: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: Initialize program
211: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212: PetscInitialize (&argc,&argv,(char*)0,help);if (ierr) return ierr;
213: MPI_Comm_size (PETSC_COMM_WORLD ,&size);
214: if (size > 1) SETERRQ (PETSC_COMM_WORLD ,PETSC_ERR_SUP,"Only for sequential runs" );
216: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
217: Create necessary matrix and vectors
218: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
219: MatCreate (PETSC_COMM_WORLD ,&A);
220: MatSetSizes (A,n,n,PETSC_DETERMINE ,PETSC_DETERMINE );
221: MatSetType (A,MATDENSE );
222: MatSetFromOptions (A);
223: MatSetUp (A);
225: MatCreateVecs (A,&U,NULL);
227: MatCreate (PETSC_COMM_WORLD ,&Jacp);
228: MatSetSizes (Jacp,PETSC_DECIDE ,PETSC_DECIDE ,2,1);
229: MatSetFromOptions (Jacp);
230: MatSetUp (Jacp);
232: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
233: Set runtime options
234: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: PetscOptionsBegin (PETSC_COMM_WORLD ,NULL,"Swing equation options" ,"" );
236: {
237: ctx.beta = 2;
238: ctx.c = 10000.0;
239: ctx.u_s = 1.0;
240: ctx.omega_s = 1.0;
241: ctx.omega_b = 120.0*PETSC_PI;
242: ctx.H = 5.0;
243: PetscOptionsScalar ("-Inertia" ,"" ,"" ,ctx.H,&ctx.H,NULL);
244: ctx.D = 5.0;
245: PetscOptionsScalar ("-D" ,"" ,"" ,ctx.D,&ctx.D,NULL);
246: ctx.E = 1.1378;
247: ctx.V = 1.0;
248: ctx.X = 0.545;
249: ctx.Pmax = ctx.E*ctx.V/ctx.X;;
250: PetscOptionsScalar ("-Pmax" ,"" ,"" ,ctx.Pmax,&ctx.Pmax,NULL);
251: ctx.Pm = 1.1;
252: PetscOptionsScalar ("-Pm" ,"" ,"" ,ctx.Pm,&ctx.Pm,NULL);
253: ctx.tf = 0.1;
254: ctx.tcl = 0.2;
255: PetscOptionsReal ("-tf" ,"Time to start fault" ,"" ,ctx.tf,&ctx.tf,NULL);
256: PetscOptionsReal ("-tcl" ,"Time to end fault" ,"" ,ctx.tcl,&ctx.tcl,NULL);
257: PetscOptionsBool ("-ensemble" ,"Run ensemble of different initial conditions" ,"" ,ensemble,&ensemble,NULL);
258: if (ensemble) {
259: ctx.tf = -1;
260: ctx.tcl = -1;
261: }
263: VecGetArray (U,&u);
264: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
265: u[1] = 1.0;
266: PetscOptionsRealArray ("-u" ,"Initial solution" ,"" ,u,&n,&flg1);
267: n = 2;
268: PetscOptionsRealArray ("-du" ,"Perturbation in initial solution" ,"" ,du,&n,&flg2);
269: u[0] += du[0];
270: u[1] += du[1];
271: VecRestoreArray (U,&u);
272: if (flg1 || flg2) {
273: ctx.tf = -1;
274: ctx.tcl = -1;
275: }
276: }
277: PetscOptionsEnd ();
279: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: Create timestepping solver context
281: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
282: TSCreate (PETSC_COMM_WORLD ,&ts);
283: TSSetProblemType (ts,TS_NONLINEAR );
284: TSSetType (ts,TSRK );
285: TSSetRHSFunction (ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
286: TSSetRHSJacobian (ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);
288: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289: Set initial conditions
290: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291: TSSetSolution (ts,U);
293: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
294: Save trajectory of solution so that TSAdjointSolve () may be used
295: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
296: TSSetSaveTrajectory (ts);
298: MatCreateVecs (A,&lambda[0],NULL);
299: /* Set initial conditions for the adjoint integration */
300: VecGetArray (lambda[0],&y_ptr);
301: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
302: VecRestoreArray (lambda[0],&y_ptr);
304: MatCreateVecs (Jacp,&mu[0],NULL);
305: VecGetArray (mu[0],&x_ptr);
306: x_ptr[0] = -1.0;
307: VecRestoreArray (mu[0],&x_ptr);
308: TSSetCostGradients (ts,1,lambda,mu);
309: TSSetCostIntegrand (ts,1,NULL,(PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec ,void*))CostIntegrand,
310: (PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec *,void*))DRDYFunction,
311: (PetscErrorCode (*)(TS ,PetscReal ,Vec ,Vec *,void*))DRDPFunction,PETSC_TRUE ,&ctx);
313: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
314: Set solver options
315: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
316: TSSetMaxTime (ts,10.0);
317: TSSetExactFinalTime (ts,TS_EXACTFINALTIME_MATCHSTEP );
318: TSSetTimeStep (ts,.01);
319: TSSetFromOptions (ts);
321: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
322: Solve nonlinear system
323: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
324: if (ensemble) {
325: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
326: VecGetArray (U,&u);
327: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
328: u[1] = ctx.omega_s;
329: u[0] += du[0];
330: u[1] += du[1];
331: VecRestoreArray (U,&u);
332: TSSetTimeStep (ts,.01);
333: TSSolve (ts,U);
334: }
335: } else {
336: TSSolve (ts,U);
337: }
338: VecView (U,PETSC_VIEWER_STDOUT_WORLD );
339: TSGetSolveTime (ts,&ftime);
340: TSGetStepNumber (ts,&steps);
342: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
343: Adjoint model starts here
344: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
345: /* Set initial conditions for the adjoint integration */
346: VecGetArray (lambda[0],&y_ptr);
347: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
348: VecRestoreArray (lambda[0],&y_ptr);
350: VecGetArray (mu[0],&x_ptr);
351: x_ptr[0] = -1.0;
352: VecRestoreArray (mu[0],&x_ptr);
354: /* Set RHS JacobianP */
355: TSSetRHSJacobianP (ts,Jacp,RHSJacobianP,&ctx);
357: TSAdjointSolve (ts);
359: PetscPrintf (PETSC_COMM_WORLD ,"\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n" );
360: VecView (lambda[0],PETSC_VIEWER_STDOUT_WORLD );
361: VecView (mu[0],PETSC_VIEWER_STDOUT_WORLD );
362: TSGetCostIntegral (ts,&q);
363: VecView (q,PETSC_VIEWER_STDOUT_WORLD );
364: VecGetArray (q,&x_ptr);
365: PetscPrintf (PETSC_COMM_WORLD ,"\n cost function=%g\n" ,(double)(x_ptr[0]-ctx.Pm));
366: VecRestoreArray (q,&x_ptr);
368: ComputeSensiP(lambda[0],mu[0],&ctx);
370: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
371: Free work space. All PETSc objects should be destroyed when they are no longer needed.
372: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
373: MatDestroy (&A);
374: MatDestroy (&Jacp);
375: VecDestroy (&U);
376: VecDestroy (&lambda[0]);
377: VecDestroy (&mu[0]);
378: TSDestroy (&ts);
379: PetscFinalize ();
380: return ierr;
381: }
384: /*TEST
386: build:
387: requires: !complex
389: test:
390: args: -viewer_binary_skip_info -ts_adapt_type none
392: TEST*/