Actual source code: ex9adj.c
petsc-3.11.1 2019-04-12
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*/