Actual source code: ex3sa.c
petsc-3.11.3 2019-06-26
2: static char help[] = "Adjoint and tangent linear sensitivity analysis of the 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}
13: /*
14: This code demonstrate the sensitivity analysis interface to a system of ordinary differential equations with discontinuities.
15: It computes the sensitivities of an integral cost function
16: \int c*max(0,\theta(t)-u_s)^beta dt
17: w.r.t. initial conditions and the parameter P_m.
18: Backward Euler method is used for time integration.
19: The discontinuities are detected with TSEvent.
20: */
22: #include <petscts.h>
24: typedef struct {
25: PetscScalar H,D,omega_b,omega_s,Pmax,Pmax_ini,Pm,E,V,X,u_s,c;
26: PetscInt beta;
27: PetscReal tf,tcl;
28: } AppCtx;
30: typedef enum {SA_ADJ, SA_TLM} SAMethod;
31: static const char *const SAMethods[] = {"ADJ","TLM","SAMethod","SA_",0};
33: /* Event check */
34: PetscErrorCode EventFunction(TS ts,PetscReal t,Vec X,PetscScalar *fvalue,void *ctx)
35: {
36: AppCtx *user=(AppCtx*)ctx;
39: /* Event for fault-on time */
40: fvalue[0] = t - user->tf;
41: /* Event for fault-off time */
42: fvalue[1] = t - user->tcl;
44: return(0);
45: }
47: PetscErrorCode PostEventFunction(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec X,PetscBool forwardsolve,void* ctx)
48: {
49: AppCtx *user=(AppCtx*)ctx;
53: if (event_list[0] == 0) {
54: if (forwardsolve) user->Pmax = 0.0; /* Apply disturbance - this is done by setting Pmax = 0 */
55: else user->Pmax = user->Pmax_ini; /* Going backward, reversal of event */
56: } else if(event_list[0] == 1) {
57: if (forwardsolve) user->Pmax = user->Pmax_ini; /* Remove the fault - this is done by setting Pmax = Pmax_ini */
58: else user->Pmax = 0.0; /* Going backward, reversal of event */
59: }
60: return(0);
61: }
63: PetscErrorCode PostStepFunction(TS ts)
64: {
65: PetscErrorCode ierr;
66: Vec U;
67: PetscReal t;
68: const PetscScalar *u;
71: TSGetTime(ts,&t);
72: TSGetSolution(ts,&U);
73: VecGetArrayRead(U,&u);
74: PetscPrintf(PETSC_COMM_SELF,"delta(%3.2f) = %8.7f\n",(double)t,(double)u[0]);
75: VecRestoreArrayRead(U,&u);
77: return(0);
78: }
80: /*
81: Defines the ODE passed to the ODE solver
82: */
83: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx *ctx)
84: {
85: PetscErrorCode ierr;
86: PetscScalar *f,Pmax;
87: const PetscScalar *u;
90: /* The next three lines allow us to access the entries of the vectors directly */
91: VecGetArrayRead(U,&u);
92: VecGetArray(F,&f);
93: Pmax = ctx->Pmax;
94: f[0] = ctx->omega_b*(u[1] - ctx->omega_s);
95: f[1] = ctx->omega_s/(2.0*ctx->H)*(ctx->Pm - Pmax*PetscSinScalar(u[0]) - ctx->D*(u[1] - ctx->omega_s));
97: VecRestoreArrayRead(U,&u);
98: VecRestoreArray(F,&f);
99: return(0);
100: }
102: /*
103: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetRHSJacobian() for the meaning of a and the Jacobian.
104: */
105: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx *ctx)
106: {
107: PetscErrorCode ierr;
108: PetscInt rowcol[] = {0,1};
109: PetscScalar J[2][2],Pmax;
110: const PetscScalar *u;
113: VecGetArrayRead(U,&u);
114: Pmax = ctx->Pmax;
115: J[0][0] = 0;
116: J[0][1] = ctx->omega_b;
117: J[1][0] = -ctx->omega_s/(2.0*ctx->H)*Pmax*PetscCosScalar(u[0]);
118: J[1][1] = -ctx->omega_s/(2.0*ctx->H)*ctx->D;
119: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
120: VecRestoreArrayRead(U,&u);
122: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
123: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
124: if (A != B) {
125: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
127: }
128: return(0);
129: }
131: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx0)
132: {
134: PetscInt row[] = {0,1},col[] = {0};
135: PetscScalar *x,J[2][1];
136: AppCtx *ctx = (AppCtx*)ctx0;
139: VecGetArray(X,&x);
140: J[0][0] = 0;
141: J[1][0] = ctx->omega_s/(2.0*ctx->H);
142: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
144: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
145: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
146: return(0);
147: }
149: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx *ctx)
150: {
151: PetscErrorCode ierr;
152: PetscScalar *r;
153: const PetscScalar *u;
156: VecGetArrayRead(U,&u);
157: VecGetArray(R,&r);
158: r[0] = ctx->c*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta);
159: VecRestoreArray(R,&r);
160: VecRestoreArrayRead(U,&u);
161: return(0);
162: }
164: static PetscErrorCode DRDYFunction(TS ts,PetscReal t,Vec U,Vec *drdy,AppCtx *ctx)
165: {
166: PetscErrorCode ierr;
167: PetscScalar *ry;
168: const PetscScalar *u;
171: VecGetArrayRead(U,&u);
172: VecGetArray(drdy[0],&ry);
173: ry[0] = ctx->c*ctx->beta*PetscPowScalarInt(PetscMax(0., u[0]-ctx->u_s),ctx->beta-1);
174: VecRestoreArray(drdy[0],&ry);
175: VecRestoreArrayRead(U,&u);
176: return(0);
177: }
179: static PetscErrorCode DRDPFunction(TS ts,PetscReal t,Vec U,Vec *drdp,AppCtx *ctx)
180: {
181: PetscErrorCode ierr;
182: PetscScalar *rp;
183: const PetscScalar *u;
186: VecGetArrayRead(U,&u);
187: VecGetArray(drdp[0],&rp);
188: rp[0] = 0.;
189: VecRestoreArray(drdp[0],&rp);
190: VecGetArrayRead(U,&u);
191: return(0);
192: }
194: PetscErrorCode ComputeSensiP(Vec lambda,Vec mu,PetscScalar *val,AppCtx *ctx)
195: {
196: PetscErrorCode ierr;
197: const PetscScalar *x,*y;
200: VecGetArrayRead(lambda,&x);
201: VecGetArrayRead(mu,&y);
202: val[0] = 1./PetscSqrtScalar(1.-(ctx->Pm/ctx->Pmax)*(ctx->Pm/ctx->Pmax))/ctx->Pmax*x[0]+y[0];
203: VecRestoreArrayRead(lambda,&x);
204: VecRestoreArrayRead(mu,&y);
205: return(0);
206: }
208: int main(int argc,char **argv)
209: {
210: TS ts; /* ODE integrator */
211: Vec U; /* solution will be stored here */
212: Mat A; /* Jacobian matrix */
213: Mat Jacp; /* JacobianP matrix */
215: PetscMPIInt size;
216: PetscInt n = 2;
217: AppCtx ctx;
218: PetscScalar *u;
219: PetscReal du[2] = {0.0,0.0};
220: PetscBool ensemble = PETSC_FALSE,flg1,flg2;
221: PetscReal ftime;
222: PetscInt steps;
223: PetscScalar *x_ptr,*y_ptr,*s_ptr;
224: Vec lambda[1],q,mu[1];
225: PetscInt direction[2];
226: PetscBool terminate[2];
227: Vec qgrad[1]; /* Forward sesivitiy */
228: Mat sp; /* Forward sensitivity matrix */
229: SAMethod sa;
231: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
232: Initialize program
233: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
235: MPI_Comm_size(PETSC_COMM_WORLD,&size);
236: if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
238: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239: Create necessary matrix and vectors
240: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
241: MatCreate(PETSC_COMM_WORLD,&A);
242: MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);
243: MatSetType(A,MATDENSE);
244: MatSetFromOptions(A);
245: MatSetUp(A);
247: MatCreateVecs(A,&U,NULL);
249: MatCreate(PETSC_COMM_WORLD,&Jacp);
250: MatSetSizes(Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
251: MatSetFromOptions(Jacp);
252: MatSetUp(Jacp);
254: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: Set runtime options
256: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
257: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");
258: {
259: ctx.beta = 2;
260: ctx.c = 10000.0;
261: ctx.u_s = 1.0;
262: ctx.omega_s = 1.0;
263: ctx.omega_b = 120.0*PETSC_PI;
264: ctx.H = 5.0;
265: PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);
266: ctx.D = 5.0;
267: PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);
268: ctx.E = 1.1378;
269: ctx.V = 1.0;
270: ctx.X = 0.545;
271: ctx.Pmax = ctx.E*ctx.V/ctx.X;
272: ctx.Pmax_ini = ctx.Pmax;
273: PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);
274: ctx.Pm = 1.1;
275: PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);
276: ctx.tf = 0.1;
277: ctx.tcl = 0.2;
278: PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);
279: PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);
280: PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL);
281: if (ensemble) {
282: ctx.tf = -1;
283: ctx.tcl = -1;
284: }
286: VecGetArray(U,&u);
287: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
288: u[1] = 1.0;
289: PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1);
290: n = 2;
291: PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2);
292: u[0] += du[0];
293: u[1] += du[1];
294: VecRestoreArray(U,&u);
295: if (flg1 || flg2) {
296: ctx.tf = -1;
297: ctx.tcl = -1;
298: }
299: sa = SA_ADJ;
300: PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)sa,(PetscEnum*)&sa,NULL);
301: }
302: PetscOptionsEnd();
304: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
305: Create timestepping solver context
306: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
307: TSCreate(PETSC_COMM_WORLD,&ts);
308: TSSetProblemType(ts,TS_NONLINEAR);
309: TSSetType(ts,TSBEULER);
310: TSSetRHSFunction(ts,NULL,(TSRHSFunction) RHSFunction,&ctx);
311: TSSetRHSJacobian(ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);
313: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
314: Set initial conditions
315: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
316: TSSetSolution(ts,U);
318: /* Set RHS JacobianP */
319: TSSetRHSJacobianP(ts,Jacp,RHSJacobianP,&ctx);
320: if (sa == SA_ADJ) {
321: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
322: Save trajectory of solution so that TSAdjointSolve() may be used
323: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
324: TSSetSaveTrajectory(ts);
326: MatCreateVecs(A,&lambda[0],NULL);
327: MatCreateVecs(Jacp,&mu[0],NULL);
328: TSSetCostGradients(ts,1,lambda,mu);
329: TSSetCostIntegrand(ts,1,NULL,(PetscErrorCode (*)(TS,PetscReal,Vec,Vec,void*))CostIntegrand,
330: (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDYFunction,
331: (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDPFunction,PETSC_TRUE,&ctx);
332: }
334: if (sa == SA_TLM) {
335: PetscScalar val[2];
336: PetscInt row[]={0,1},col[]={0};
338: VecCreate(PETSC_COMM_WORLD,&qgrad[0]);
339: VecSetSizes(qgrad[0],PETSC_DECIDE,1);
340: VecSetFromOptions(qgrad[0]);
342: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp);
343: TSForwardSetSensitivities(ts,1,sp);
344: TSForwardSetIntegralGradients(ts,1,qgrad);
345: TSSetCostIntegrand(ts,1,NULL,(PetscErrorCode (*)(TS,PetscReal,Vec,Vec,void*))CostIntegrand,
346: (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDYFunction,
347: (PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDPFunction,PETSC_TRUE,&ctx);
348: val[0] = 1./PetscSqrtScalar(1.-(ctx.Pm/ctx.Pmax)*(ctx.Pm/ctx.Pmax))/ctx.Pmax;
349: val[1] = 0.0;
350: MatSetValues(sp,2,row,1,col,val,INSERT_VALUES);
351: MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY);
352: MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY);
353: }
355: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356: Set solver options
357: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
358: TSSetMaxTime(ts,1.0);
359: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
360: TSSetTimeStep(ts,0.03125);
361: TSSetFromOptions(ts);
363: direction[0] = direction[1] = 1;
364: terminate[0] = terminate[1] = PETSC_FALSE;
366: TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&ctx);
368: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
369: Solve nonlinear system
370: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
371: if (ensemble) {
372: for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
373: VecGetArray(U,&u);
374: u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
375: u[1] = ctx.omega_s;
376: u[0] += du[0];
377: u[1] += du[1];
378: VecRestoreArray(U,&u);
379: TSSetTimeStep(ts,0.03125);
380: TSSolve(ts,U);
381: }
382: } else {
383: TSSolve(ts,U);
384: }
385: TSGetSolveTime(ts,&ftime);
386: TSGetStepNumber(ts,&steps);
388: if (sa == SA_ADJ) {
389: PetscScalar grad[1];
390: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
391: Adjoint model starts here
392: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
393: /* Set initial conditions for the adjoint integration */
394: VecGetArray(lambda[0],&y_ptr);
395: y_ptr[0] = 0.0; y_ptr[1] = 0.0;
396: VecRestoreArray(lambda[0],&y_ptr);
398: VecGetArray(mu[0],&x_ptr);
399: x_ptr[0] = 0.0;
400: VecRestoreArray(mu[0],&x_ptr);
402: TSAdjointSolve(ts);
404: PetscPrintf(PETSC_COMM_WORLD,"\n lambda: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n");
405: VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD);
406: PetscPrintf(PETSC_COMM_WORLD,"\n mu: d[Psi(tf)]/d[pm] d[Psi(tf)]/d[pm]\n");
407: VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD);
408: TSGetCostIntegral(ts,&q);
409: VecGetArray(q,&x_ptr);
410: PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm));
411: VecRestoreArray(q,&x_ptr);
412: ComputeSensiP(lambda[0],mu[0],grad,&ctx);
413: PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)grad[0]);
414: VecDestroy(&lambda[0]);
415: VecDestroy(&mu[0]);
416: }
418: if (sa == SA_TLM) {
419: PetscPrintf(PETSC_COMM_WORLD,"\n trajectory sensitivity: d[Psi(tf)]/d[pm] d[Psi(tf)]/d[pm]\n");
420: MatView(sp,PETSC_VIEWER_STDOUT_WORLD);
421: TSGetCostIntegral(ts,&q);
422: VecGetArray(q,&s_ptr);
423: PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(s_ptr[0]-ctx.Pm));
424: VecRestoreArray(q,&s_ptr);
425: VecGetArray(qgrad[0],&s_ptr);
426: PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)s_ptr[0]);
427: VecRestoreArray(qgrad[0],&s_ptr);
428: VecDestroy(&qgrad[0]);
429: MatDestroy(&sp);
430: }
431: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432: Free work space. All PETSc objects should be destroyed when they are no longer needed.
433: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
434: MatDestroy(&A);
435: MatDestroy(&Jacp);
436: VecDestroy(&U);
437: TSDestroy(&ts);
438: PetscFinalize();
439: return ierr;
440: }
443: /*TEST
445: build:
446: requires: !complex !single
448: test:
449: args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu
451: test:
452: suffix: 2
453: args: -sa_method tlm -ts_type cn -pc_type lu
455: test:
456: suffix: 3
457: args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp
459: TEST*/