Actual source code: ex9adj.c
petsc-3.12.0 2019-09-29
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: return (0);
57: }
59: /*
60: Defines the ODE passed to the ODE solver
61: */
62: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
63: {
64: PetscErrorCode ierr;
65: PetscScalar *f,Pmax;
66: const PetscScalar *u;
69: /* The next three lines allow us to access the entries of the vectors directly */
70: VecGetArrayRead (U,&u);
71: VecGetArray (F,&f);
72: 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 */
73: else Pmax = ctx->Pmax;
75: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
76: f[1] = (-Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s) + ctx->Pm)*ctx->omega_s/(2.0*ctx->H);
78: VecRestoreArrayRead (U,&u);
79: VecRestoreArray (F,&f);
80: return (0);
81: }
83: /*
84: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian () for the meaning of a and the Jacobian.
85: */
86: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
87: {
88: PetscErrorCode ierr;
89: PetscInt rowcol[] = {0,1};
90: PetscScalar J[2][2],Pmax;
91: const PetscScalar *u;
94: VecGetArrayRead (U,&u);
95: 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 */
96: else Pmax = ctx->Pmax;
98: J[0][0] = 0; J[0][1] = ctx->omega_b;
99: 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);
101: MatSetValues (A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES );
102: VecRestoreArrayRead (U,&u);
104: MatAssemblyBegin (A,MAT_FINAL_ASSEMBLY );
105: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY );
106: if (A != B) {
107: MatAssemblyBegin (B,MAT_FINAL_ASSEMBLY );
108: MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY );
109: }
110: return (0);
111: }
113: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
114: {
116: PetscInt row[] = {0,1},col[]={0};
117: PetscScalar J[2][1];
118: AppCtx *ctx=(AppCtx*)ctx0;
121: J[0][0] = 0;
122: J[1][0] = ctx->omega_s/(2.0*ctx->H);
123: MatSetValues (A,2,row,1,col,&J[0][0],INSERT_VALUES );
124: MatAssemblyBegin (A,MAT_FINAL_ASSEMBLY );
125: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY );
126: return (0);
127: }
129: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
130: {
131: PetscErrorCode ierr;
132: PetscScalar *r;
133: const PetscScalar *u;
136: VecGetArrayRead (U,&u);
137: VecGetArray (R,&r);
138: r[0] = ctx->c*PetscPowScalarInt(PetscMax (0., u[0]-ctx->u_s),ctx->beta);
139: VecRestoreArray (R,&r);
140: VecRestoreArrayRead (U,&u);
141: return (0);
142: }
144: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx *ctx)
145: {
146: PetscErrorCode ierr;
147: PetscScalar ru[1];
148: const PetscScalar *u;
149: PetscInt row[] = {0},col[] = {0};
152: VecGetArrayRead (U,&u);
153: ru[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax (0., u[0]-ctx->u_s),ctx->beta-1);
154: VecRestoreArrayRead (U,&u);
155: MatSetValues (DRDU,1,row,1,col,ru,INSERT_VALUES );
156: MatAssemblyBegin (DRDU,MAT_FINAL_ASSEMBLY );
157: MatAssemblyEnd (DRDU,MAT_FINAL_ASSEMBLY );
158: return (0);
159: }
161: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx *ctx)
162: {
163: PetscErrorCode ierr;
166: MatZeroEntries (DRDP);
167: MatAssemblyBegin (DRDP,MAT_FINAL_ASSEMBLY );
168: MatAssemblyEnd (DRDP,MAT_FINAL_ASSEMBLY );
169: return (0);
170: }
172: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,AppCtx *ctx)
173: {
174: PetscErrorCode ierr;
175: PetscScalar sensip;
176: const PetscScalar *x,*y;
179: VecGetArrayRead (lambda,&x);
180: VecGetArrayRead (mu,&y);
181: sensip = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
182: PetscPrintf (PETSC_COMM_WORLD ,"\n sensitivity wrt parameter pm: %.7f \n" ,(double)sensip);
183: VecRestoreArrayRead (lambda,&x);
184: VecRestoreArrayRead (mu,&y);
185: return (0);
186: }
188: int main(int argc,char **argv)
189: {
190: TS ts,quadts; /* ODE integrator */
191: Vec U; /* solution will be stored here */
192: Mat A; /* Jacobian matrix */
193: Mat Jacp; /* Jacobian matrix */
194: Mat DRDU,DRDP;
196: PetscMPIInt size;
197: PetscInt n = 2;
198: AppCtx ctx;
199: PetscScalar *u;
200: PetscReal du[2] = {0.0,0.0};
201: PetscBool ensemble = PETSC_FALSE ,flg1,flg2;
202: PetscReal ftime;
203: PetscInt steps;
204: PetscScalar *x_ptr,*y_ptr;
205: Vec lambda[1],q,mu[1];
207: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208: Initialize program
209: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
210: PetscInitialize (&argc,&argv,(char*)0,help);if (ierr) return ierr;
211: MPI_Comm_size (PETSC_COMM_WORLD ,&size);
212: if (size > 1) SETERRQ (PETSC_COMM_WORLD ,PETSC_ERR_SUP,"Only for sequential runs" );
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Create necessary matrix and vectors
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: MatCreate (PETSC_COMM_WORLD ,&A);
218: MatSetSizes (A,n,n,PETSC_DETERMINE ,PETSC_DETERMINE );
219: MatSetType (A,MATDENSE );
220: MatSetFromOptions (A);
221: MatSetUp (A);
223: MatCreateVecs (A,&U,NULL);
225: MatCreate (PETSC_COMM_WORLD ,&Jacp);
226: MatSetSizes (Jacp,PETSC_DECIDE ,PETSC_DECIDE ,2,1);
227: MatSetFromOptions (Jacp);
228: MatSetUp (Jacp);
230: MatCreateDense (PETSC_COMM_WORLD ,PETSC_DECIDE ,PETSC_DECIDE ,1,1,NULL,&DRDP);
231: MatSetUp (DRDP);
232: MatCreateDense (PETSC_COMM_WORLD ,PETSC_DECIDE ,PETSC_DECIDE ,1,1,NULL,&DRDU);
233: MatSetUp (DRDU);
235: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236: Set runtime options
237: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238: PetscOptionsBegin (PETSC_COMM_WORLD ,NULL,"Swing equation options" ,"" );
239: {
240: ctx.beta = 2;
241: ctx.c = 10000.0;
242: ctx.u_s = 1.0;
243: ctx.omega_s = 1.0;
244: ctx.omega_b = 120.0*PETSC_PI;
245: ctx.H = 5.0;
246: PetscOptionsScalar ("-Inertia" ,"" ,"" ,ctx.H,&ctx.H,NULL);
247: ctx.D = 5.0;
248: PetscOptionsScalar ("-D" ,"" ,"" ,ctx.D,&ctx.D,NULL);
249: ctx.E = 1.1378;
250: ctx.V = 1.0;
251: ctx.X = 0.545;
252: ctx.Pmax = ctx.E*ctx.V/ctx.X;;
253: PetscOptionsScalar ("-Pmax" ,"" ,"" ,ctx.Pmax,&ctx.Pmax,NULL);
254: ctx.Pm = 1.1;
255: PetscOptionsScalar ("-Pm" ,"" ,"" ,ctx.Pm,&ctx.Pm,NULL);
256: ctx.tf = 0.1;
257: ctx.tcl = 0.2;
258: PetscOptionsReal ("-tf" ,"Time to start fault" ,"" ,ctx.tf,&ctx.tf,NULL);
259: PetscOptionsReal ("-tcl" ,"Time to end fault" ,"" ,ctx.tcl,&ctx.tcl,NULL);
260: PetscOptionsBool ("-ensemble" ,"Run ensemble of different initial conditions" ,"" ,ensemble,&ensemble,NULL);
261: if (ensemble) {
262: ctx.tf = -1;
263: ctx.tcl = -1;
264: }
266: VecGetArray (U,&u);
267: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
268: u[1] = 1.0;
269: PetscOptionsRealArray ("-u" ,"Initial solution" ,"" ,u,&n,&flg1);
270: n = 2;
271: PetscOptionsRealArray ("-du" ,"Perturbation in initial solution" ,"" ,du,&n,&flg2);
272: u[0] += du[0];
273: u[1] += du[1];
274: VecRestoreArray (U,&u);
275: if (flg1 || flg2) {
276: ctx.tf = -1;
277: ctx.tcl = -1;
278: }
279: }
280: PetscOptionsEnd ();
282: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
283: Create timestepping solver context
284: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
285: TSCreate (PETSC_COMM_WORLD ,&ts);
286: TSSetProblemType (ts,TS_NONLINEAR );
287: TSSetEquationType (ts,TS_EQ_ODE_EXPLICIT ); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
288: TSSetType (ts,TSRK );
289: TSSetRHSFunction (ts,NULL,(TSRHSFunction)RHSFunction,&ctx);
290: TSSetRHSJacobian (ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);
291: TSCreateQuadratureTS (ts,PETSC_TRUE ,&quadts);
292: TSSetRHSFunction (quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);
293: TSSetRHSJacobian (quadts,DRDU,DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);
294: TSSetRHSJacobianP (quadts,DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&ctx);
296: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
297: Set initial conditions
298: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
299: TSSetSolution (ts,U);
301: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
302: Save trajectory of solution so that TSAdjointSolve () may be used
303: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
304: TSSetSaveTrajectory (ts);
306: MatCreateVecs (A,&lambda[0],NULL);
307: /* Set initial conditions for the adjoint integration */
308: VecGetArray (lambda[0],&y_ptr);
309: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
310: VecRestoreArray (lambda[0],&y_ptr);
312: MatCreateVecs (Jacp,&mu[0],NULL);
313: VecGetArray (mu[0],&x_ptr);
314: x_ptr[0] = -1.0;
315: VecRestoreArray (mu[0],&x_ptr);
316: TSSetCostGradients (ts,1,lambda,mu);
318: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
319: Set solver options
320: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
321: TSSetMaxTime (ts,10.0);
322: TSSetExactFinalTime (ts,TS_EXACTFINALTIME_MATCHSTEP );
323: TSSetTimeStep (ts,.01);
324: TSSetFromOptions (ts);
326: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
327: Solve nonlinear system
328: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
329: if (ensemble) {
330: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
331: VecGetArray (U,&u);
332: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
333: u[1] = ctx.omega_s;
334: u[0] += du[0];
335: u[1] += du[1];
336: VecRestoreArray (U,&u);
337: TSSetTimeStep (ts,.01);
338: TSSolve (ts,U);
339: }
340: } else {
341: TSSolve (ts,U);
342: }
343: VecView (U,PETSC_VIEWER_STDOUT_WORLD );
344: TSGetSolveTime (ts,&ftime);
345: TSGetStepNumber (ts,&steps);
347: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: Adjoint model starts here
349: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350: /* Set initial conditions for the adjoint integration */
351: VecGetArray (lambda[0],&y_ptr);
352: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
353: VecRestoreArray (lambda[0],&y_ptr);
355: VecGetArray (mu[0],&x_ptr);
356: x_ptr[0] = -1.0;
357: VecRestoreArray (mu[0],&x_ptr);
359: /* Set RHS JacobianP */
360: TSSetRHSJacobianP (ts,Jacp,RHSJacobianP,&ctx);
362: TSAdjointSolve (ts);
364: PetscPrintf (PETSC_COMM_WORLD ,"\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n" );
365: VecView (lambda[0],PETSC_VIEWER_STDOUT_WORLD );
366: VecView (mu[0],PETSC_VIEWER_STDOUT_WORLD );
367: TSGetCostIntegral (ts,&q);
368: VecView (q,PETSC_VIEWER_STDOUT_WORLD );
369: VecGetArray (q,&x_ptr);
370: PetscPrintf (PETSC_COMM_WORLD ,"\n cost function=%g\n" ,(double)(x_ptr[0]-ctx.Pm));
371: VecRestoreArray (q,&x_ptr);
373: ComputeSensiP(lambda[0],mu[0],&ctx);
375: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
376: Free work space. All PETSc objects should be destroyed when they are no longer needed.
377: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
378: MatDestroy (&A);
379: MatDestroy (&Jacp);
380: MatDestroy (&DRDU);
381: MatDestroy (&DRDP);
382: VecDestroy (&U);
383: VecDestroy (&lambda[0]);
384: VecDestroy (&mu[0]);
385: TSDestroy (&ts);
386: PetscFinalize ();
387: return ierr;
388: }
391: /*TEST
393: build:
394: requires: !complex
396: test:
397: args: -viewer_binary_skip_info -ts_adapt_type none
399: TEST*/