1
2 static char help[] = "Adjoint and tangent linear sensitivity analysis of the basic equation for generator stability analysis.\n";
3
4 /*F
5
6 \begin{eqnarray}
7 \frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
8 \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
9 \end{eqnarray}
10
11 F*/
12
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 */
21
22 #include <petscts.h>
23 #include "ex3.h"
24
main(int argc,char ** argv)25 int main(int argc,char **argv)
26 {
27 TS ts,quadts; /* ODE integrator */
28 Vec U; /* solution will be stored here */
29 PetscErrorCode ierr;
30 PetscMPIInt size;
31 PetscInt n = 2;
32 AppCtx ctx;
33 PetscScalar *u;
34 PetscReal du[2] = {0.0,0.0};
35 PetscBool ensemble = PETSC_FALSE,flg1,flg2;
36 PetscReal ftime;
37 PetscInt steps;
38 PetscScalar *x_ptr,*y_ptr,*s_ptr;
39 Vec lambda[1],q,mu[1];
40 PetscInt direction[2];
41 PetscBool terminate[2];
42 Mat qgrad;
43 Mat sp; /* Forward sensitivity matrix */
44 SAMethod sa;
45
46 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47 Initialize program
48 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
49 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
50 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
51 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs");
52
53 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54 Create necessary matrix and vectors
55 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
56 ierr = MatCreate(PETSC_COMM_WORLD,&ctx.Jac);CHKERRQ(ierr);
57 ierr = MatSetSizes(ctx.Jac,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
58 ierr = MatSetType(ctx.Jac,MATDENSE);CHKERRQ(ierr);
59 ierr = MatSetFromOptions(ctx.Jac);CHKERRQ(ierr);
60 ierr = MatSetUp(ctx.Jac);CHKERRQ(ierr);
61 ierr = MatCreateVecs(ctx.Jac,&U,NULL);CHKERRQ(ierr);
62 ierr = MatCreate(PETSC_COMM_WORLD,&ctx.Jacp);CHKERRQ(ierr);
63 ierr = MatSetSizes(ctx.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);CHKERRQ(ierr);
64 ierr = MatSetFromOptions(ctx.Jacp);CHKERRQ(ierr);
65 ierr = MatSetUp(ctx.Jacp);CHKERRQ(ierr);
66 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&ctx.DRDP);CHKERRQ(ierr);
67 ierr = MatSetUp(ctx.DRDP);CHKERRQ(ierr);
68 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&ctx.DRDU);CHKERRQ(ierr);
69 ierr = MatSetUp(ctx.DRDU);CHKERRQ(ierr);
70
71 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72 Set runtime options
73 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Swing equation options","");CHKERRQ(ierr);
75 {
76 ctx.beta = 2;
77 ctx.c = 10000.0;
78 ctx.u_s = 1.0;
79 ctx.omega_s = 1.0;
80 ctx.omega_b = 120.0*PETSC_PI;
81 ctx.H = 5.0;
82 ierr = PetscOptionsScalar("-Inertia","","",ctx.H,&ctx.H,NULL);CHKERRQ(ierr);
83 ctx.D = 5.0;
84 ierr = PetscOptionsScalar("-D","","",ctx.D,&ctx.D,NULL);CHKERRQ(ierr);
85 ctx.E = 1.1378;
86 ctx.V = 1.0;
87 ctx.X = 0.545;
88 ctx.Pmax = ctx.E*ctx.V/ctx.X;
89 ctx.Pmax_ini = ctx.Pmax;
90 ierr = PetscOptionsScalar("-Pmax","","",ctx.Pmax,&ctx.Pmax,NULL);CHKERRQ(ierr);
91 ctx.Pm = 1.1;
92 ierr = PetscOptionsScalar("-Pm","","",ctx.Pm,&ctx.Pm,NULL);CHKERRQ(ierr);
93 ctx.tf = 0.1;
94 ctx.tcl = 0.2;
95 ierr = PetscOptionsReal("-tf","Time to start fault","",ctx.tf,&ctx.tf,NULL);CHKERRQ(ierr);
96 ierr = PetscOptionsReal("-tcl","Time to end fault","",ctx.tcl,&ctx.tcl,NULL);CHKERRQ(ierr);
97 ierr = PetscOptionsBool("-ensemble","Run ensemble of different initial conditions","",ensemble,&ensemble,NULL);CHKERRQ(ierr);
98 if (ensemble) {
99 ctx.tf = -1;
100 ctx.tcl = -1;
101 }
102
103 ierr = VecGetArray(U,&u);CHKERRQ(ierr);
104 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
105 u[1] = 1.0;
106 ierr = PetscOptionsRealArray("-u","Initial solution","",u,&n,&flg1);CHKERRQ(ierr);
107 n = 2;
108 ierr = PetscOptionsRealArray("-du","Perturbation in initial solution","",du,&n,&flg2);CHKERRQ(ierr);
109 u[0] += du[0];
110 u[1] += du[1];
111 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
112 if (flg1 || flg2) {
113 ctx.tf = -1;
114 ctx.tcl = -1;
115 }
116 sa = SA_ADJ;
117 ierr = PetscOptionsEnum("-sa_method","Sensitivity analysis method (adj or tlm)","",SAMethods,(PetscEnum)sa,(PetscEnum*)&sa,NULL);CHKERRQ(ierr);
118 }
119 ierr = PetscOptionsEnd();CHKERRQ(ierr);
120
121 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122 Create timestepping solver context
123 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
124 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
125 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
126 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
127 ierr = TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);CHKERRQ(ierr);
128 ierr = TSSetRHSJacobian(ts,ctx.Jac,ctx.Jac,(TSRHSJacobian)RHSJacobian,&ctx);CHKERRQ(ierr);
129
130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131 Set initial conditions
132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133 ierr = TSSetSolution(ts,U);CHKERRQ(ierr);
134
135 /* Set RHS JacobianP */
136 ierr = TSSetRHSJacobianP(ts,ctx.Jacp,RHSJacobianP,&ctx);CHKERRQ(ierr);
137
138 ierr = TSCreateQuadratureTS(ts,PETSC_FALSE,&quadts);CHKERRQ(ierr);
139 ierr = TSSetRHSFunction(quadts,NULL,(TSRHSFunction)CostIntegrand,&ctx);CHKERRQ(ierr);
140 ierr = TSSetRHSJacobian(quadts,ctx.DRDU,ctx.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&ctx);CHKERRQ(ierr);
141 ierr = TSSetRHSJacobianP(quadts,ctx.DRDP,DRDPJacobianTranspose,&ctx);CHKERRQ(ierr);
142 if (sa == SA_ADJ) {
143 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144 Save trajectory of solution so that TSAdjointSolve() may be used
145 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
146 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);
147 ierr = MatCreateVecs(ctx.Jac,&lambda[0],NULL);CHKERRQ(ierr);
148 ierr = MatCreateVecs(ctx.Jacp,&mu[0],NULL);CHKERRQ(ierr);
149 ierr = TSSetCostGradients(ts,1,lambda,mu);CHKERRQ(ierr);
150 }
151
152 if (sa == SA_TLM) {
153 PetscScalar val[2];
154 PetscInt row[]={0,1},col[]={0};
155
156 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,&qgrad);CHKERRQ(ierr);
157 ierr = MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&sp);CHKERRQ(ierr);
158 ierr = TSForwardSetSensitivities(ts,1,sp);CHKERRQ(ierr);
159 ierr = TSForwardSetSensitivities(quadts,1,qgrad);CHKERRQ(ierr);
160 val[0] = 1./PetscSqrtScalar(1.-(ctx.Pm/ctx.Pmax)*(ctx.Pm/ctx.Pmax))/ctx.Pmax;
161 val[1] = 0.0;
162 ierr = MatSetValues(sp,2,row,1,col,val,INSERT_VALUES);CHKERRQ(ierr);
163 ierr = MatAssemblyBegin(sp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
164 ierr = MatAssemblyEnd(sp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
165 }
166
167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168 Set solver options
169 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170 ierr = TSSetMaxTime(ts,1.0);CHKERRQ(ierr);
171 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
172 ierr = TSSetTimeStep(ts,0.03125);CHKERRQ(ierr);
173 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
174
175 direction[0] = direction[1] = 1;
176 terminate[0] = terminate[1] = PETSC_FALSE;
177
178 ierr = TSSetEventHandler(ts,2,direction,terminate,EventFunction,PostEventFunction,(void*)&ctx);CHKERRQ(ierr);
179
180 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181 Solve nonlinear system
182 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183 if (ensemble) {
184 for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
185 ierr = VecGetArray(U,&u);CHKERRQ(ierr);
186 u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
187 u[1] = ctx.omega_s;
188 u[0] += du[0];
189 u[1] += du[1];
190 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
191 ierr = TSSetTimeStep(ts,0.03125);CHKERRQ(ierr);
192 ierr = TSSolve(ts,U);CHKERRQ(ierr);
193 }
194 } else {
195 ierr = TSSolve(ts,U);CHKERRQ(ierr);
196 }
197 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
198 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
199
200 if (sa == SA_ADJ) {
201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202 Adjoint model starts here
203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204 /* Set initial conditions for the adjoint integration */
205 ierr = VecGetArray(lambda[0],&y_ptr);CHKERRQ(ierr);
206 y_ptr[0] = 0.0; y_ptr[1] = 0.0;
207 ierr = VecRestoreArray(lambda[0],&y_ptr);CHKERRQ(ierr);
208
209 ierr = VecGetArray(mu[0],&x_ptr);CHKERRQ(ierr);
210 x_ptr[0] = 0.0;
211 ierr = VecRestoreArray(mu[0],&x_ptr);CHKERRQ(ierr);
212
213 ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
214
215 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n lambda: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n");CHKERRQ(ierr);
216 ierr = VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
217 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n mu: d[Psi(tf)]/d[pm]\n");CHKERRQ(ierr);
218 ierr = VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
219 ierr = TSGetCostIntegral(ts,&q);CHKERRQ(ierr);
220 ierr = VecGetArray(q,&x_ptr);CHKERRQ(ierr);
221 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm));CHKERRQ(ierr);
222 ierr = VecRestoreArray(q,&x_ptr);CHKERRQ(ierr);
223 ierr = ComputeSensiP(lambda[0],mu[0],&ctx);CHKERRQ(ierr);
224 ierr = VecGetArray(mu[0],&x_ptr);CHKERRQ(ierr);
225 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)x_ptr[0]);CHKERRQ(ierr);
226 ierr = VecRestoreArray(mu[0],&x_ptr);CHKERRQ(ierr);
227 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr);
228 ierr = VecDestroy(&mu[0]);CHKERRQ(ierr);
229 }
230 if (sa == SA_TLM) {
231 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n trajectory sensitivity: d[phi(tf)]/d[pm] d[omega(tf)]/d[pm]\n");CHKERRQ(ierr);
232 ierr = MatView(sp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
233 ierr = TSGetCostIntegral(ts,&q);CHKERRQ(ierr);
234 ierr = VecGetArray(q,&s_ptr);CHKERRQ(ierr);
235 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(s_ptr[0]-ctx.Pm));CHKERRQ(ierr);
236 ierr = VecRestoreArray(q,&s_ptr);CHKERRQ(ierr);
237 ierr = MatDenseGetArray(qgrad,&s_ptr);CHKERRQ(ierr);
238 ierr = PetscPrintf(PETSC_COMM_WORLD,"\n gradient=%g\n",(double)s_ptr[0]);CHKERRQ(ierr);
239 ierr = MatDenseRestoreArray(qgrad,&s_ptr);CHKERRQ(ierr);
240 ierr = MatDestroy(&qgrad);CHKERRQ(ierr);
241 ierr = MatDestroy(&sp);CHKERRQ(ierr);
242 }
243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244 Free work space. All PETSc objects should be destroyed when they are no longer needed.
245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
246 ierr = MatDestroy(&ctx.Jac);CHKERRQ(ierr);
247 ierr = MatDestroy(&ctx.Jacp);CHKERRQ(ierr);
248 ierr = MatDestroy(&ctx.DRDU);CHKERRQ(ierr);
249 ierr = MatDestroy(&ctx.DRDP);CHKERRQ(ierr);
250 ierr = VecDestroy(&U);CHKERRQ(ierr);
251 ierr = TSDestroy(&ts);CHKERRQ(ierr);
252 ierr = PetscFinalize();
253 return ierr;
254 }
255
256
257 /*TEST
258
259 build:
260 requires: !complex !single
261
262 test:
263 args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu
264
265 test:
266 suffix: 2
267 args: -sa_method tlm -ts_type cn -pc_type lu
268
269 test:
270 suffix: 3
271 args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp
272
273 TEST*/
274